First-order Raman spectra of wurtzite GaN made from natural Ga and isotopically pure ${}^{15}\mathrm{N}$ and natural N (99.63% ${}^{14}\mathrm{N}$) were measured at low temperature. The Raman frequencies of the polar optical phonons [${A}_{1}(\mathrm{TO},\mathrm{LO})$, ${E}_{1}(\mathrm{TO},\mathrm{LO})$] shift according to the inverse square root of the reduced masses, as expected. The isotope shifts of the two nonpolar ${E}_{2}$ modes deviate significantly from the reduced-mass behavior as well as from the dependence expected for pure N or Ga vibrations which might be expected in view of their large frequency difference. This indicates that the ${E}_{2}$ modes involve mixed Ga and N vibrations. From fits to the experimental data with a coupled two-mode model we determine the mode eigenvectors. They confirm the results of an ab initio calculation which we have done for the zone-center vibrations of wurtzite GaN. Our calculations also predict a vanishing coupling between the two ${B}_{1}$ silent modes. We also measured a ${}^{\mathrm{nat}}{\mathrm{Ga}}^{14}{\mathrm{N}}_{0.5}^{15}{\mathrm{N}}_{0.5}$ alloy sample. Additional negative frequency shifts due to isotope disorder are observed for the ${A}_{1}(\mathrm{TO},\mathrm{LO})$, ${E}_{1}(\mathrm{TO},\mathrm{LO})$, and the high-frequency ${E}_{2}$ modes as compared to the reduced-mass or coupled-mode behavior. An estimate from second-order perturbation theory provides a quantitative explanation of this result.