Close correspondence between quantitative and moleculargenetic divergence times for Neandertals and modern humans
 *Department of Anthropology, University of California, One Shields Avenue, Davis, CA 95616;
 ^{†}Department of Human Evolution, Max Planck Institute for Evolutionary Anthropology, Deutscher Platz 6, D04103 Leipzig, Germany;
 ^{§}Department of Anthropology, University of Illinois at Urbana–Champaign, 109 Davenport Hall, 607 South Matthews Avenue, Urbana, IL 61801; and
 ^{¶}Department of Palaeontology, Natural History Museum, London SW7 5BD, United Kingdom
See allHide authors and affiliations

Edited by Erik Trinkaus, Washington University, St. Louis, MO, and approved January 17, 2008 (received for review September 27, 2007)
Abstract
Recent research has shown that genetic drift may have produced many cranial differences between Neandertals and modern humans. If this is the case, then it should be possible to estimate population genetic parameters from Neandertal and modern human cranial measurements in a manner analogous to how estimates are made from DNA sequences. Building on previous work in evolutionary quantitative genetics and on microsatellites, we present a divergence time estimator for neutrally evolving morphological measurements. We then apply this estimator to 37 standard cranial measurements collected on 2,524 modern humans from 30 globally distributed populations and 20 Neandertal specimens. We calculate that the lineages leading to Neandertals and modern humans split ≈311,000 (95% C.I.: 182,000 to 466,000) or 435,000 (95% C.I.: 308,000 to 592,000) years ago, depending on assumptions about changes in withinpopulation variation. These dates are quite similar to those recently derived from ancient Neandertal and extant human DNA sequences. Close correspondence between cranial and DNAsequence results implies that both datasets largely, although not necessarily exclusively, reflect neutral divergence, causing them to track population history or phylogeny rather than the action of diversifying natural selection. The cranial dataset covers only aspects of cranial anatomy that can be readily quantified with standard osteometric tools, so future research will be needed to determine whether these results are representative. Nonetheless, for the measurements we consider here, we find no conflict between molecules and morphology.
Lately, evidence has been accumulating for the importance of neutral evolution in producing cranial differences among human populations and between Neandertals ^{‖} and modern humans. Under neutral evolution, genetic drift provides the mechanism, and mutation provides the raw material for divergence between groups, with natural selection being relegated to a role of slowing down divergence (1–3). Lynch (4) and Relethford (5) provided some of the first evidence of the importance of neutral evolution for understanding human cranial differences. Lynch (4) found that cranial distances among human populations corresponded well with relative and absolute rates of divergence predicted by neutral evolution. Relethford (5) showed that estimates of F _{ST} (a measure of amongpopulation differentiation) from human cranial measurements were similar to those from presumably neutral genetic loci. Further work established that cranial and molecular distances among human populations tend to be significantly associated with each other (6–9), both cranial and molecular distances are correlated with geographic distances among globally distributed human populations (10, 11), cranial measurements appear to fit neutral expectations as well as microsatellites for humans (12), and statistical tests fail to detect deviations from neutrality for cranial differences between Neandertals and modern humans, even though the tests look reasonably powerful (12). Furthermore, withinpopulation cranial diversity seems to decrease with geographic distance from subSaharan Africa (13). This decrease may mirror a decrease in molecular diversity (14–17), but the strength of the cranial relationship is much weaker.
If neutral evolution is responsible for cranial differences among human populations and between Neandertals and modern humans, then it should be possible to estimate population genetic parameters from cranial measurements in a manner analogous to how estimates are made from DNA sequences. Building on previous work in evolutionary quantitative genetics (18–21) and on microsatellites (22–24), we present a divergence time estimator for neutrally evolving morphological measurements. We then apply this estimator to cranial measurements to infer when Neandertals and modern humans diverged, and we compare our results with those from ancient Neandertal and extant human DNA sequences.
Morphological Divergence Time Estimator
Divergence Time.
To estimate the divergence time of two populations (or species), we adapt the T _{D} estimator that was originally developed for microsatellites (23). Microsatellites are rapidly evolving blocks of DNA for which a simple DNA sequence is repeated multiple times, and individuals vary in their number of repeats. So, like morphological measurements (metric characteristics), microsatellites are quantitative characters that change by lengthening and shortening (24–26), making the adaptation fairly straightforward. We call the new estimator PT _{D} for phenotypic T _{D} and to distinguish it from T _{D} for microsatellites. By divergence time, we mean when the two populations last shared a randomly mating common ancestor.
We define PT _{D} as follows. For two daughter populations at mutation drift equilibrium (balance between the addition of variation by mutation and the removal of variation by genetic drift), the betweenpopulation variance for a measurement is expected to increase at the rate of 2 V _{m} per generation (20, 21), where V _{m} is the average amount of new additive genetic variance introduced by mutation per zygote per generation in each population. This result holds for many different underlying genetic models (20) as long as the divergence is by genetic drift rather than by natural selection. Let x _{1} and x _{2} be measurement means for populations nos. 1 and 2, respectively. An estimator for the time in generations when the two populations diverged, under the assumption of mutation drift equilibrium, is: Note that (x _{1}−x _{2})^{2} equals twice the betweenpopulation variance, which is the quantity most commonly used in the evolutionary quantitative genetics literature (21). Also note that the numerator in Eq. 1 is mathematically analogous to the expression for (δμ)^{2}, a genetic distance for microsatellites (22, 23), if a measurement is equated to a single microsatellite. The denominator in Eq. 1 , or the expected rate of increase, is four times the mutation parameter, V _{m}, whereas for (δμ)^{2} , it is twice the mutation parameter (22, 23). This is because of differences in the definitions of the mutation parameters. By necessity, all withinpopulation and mutational variances in the quantitative genetics literature are zygotic, whereas they are gametic for microsatellites. Gametic variances include differences both among individuals and among the paired chromosomes of a given individual; zygotic variances only include amongindividual differences. At Hardy–Weinberg equilibrium, a gametic variance should be twice a zygotic variance (27), hence the difference by a factor of two between the quantitative genetic and microsatellite formulas.
Following results for microsatellites (23), if the assumption of mutation drift equilibrium is relaxed, three additional variance terms must be added to the numerator (with modifications to account for zygotic vs. gametic variances as before). These terms are measures of the degree of departure from mutation drift equilibrium. Departures could happen with population expansions or migration between populations. Let V _{1} and V _{2} be the additive genetic variances for the measurements for populations nos. 1 and 2, respectively. Let V _{0} be the additive genetic variance in the ancestral population before subdivision into two daughter populations. Under a more general model, an estimator for the time in generations when the two populations diverged is: Let h ^{2} be the narrowsense heritability for the measurement in both populations and σ_{1} ^{2} and σ_{2} ^{2} be the phenotypic variances for populations nos. 1 and 2, respectively. Following Falconer (28), Let m be a mutation parameter and σ_{P} ^{2} be the pooled within population phenotypic variance. Following Lynch (29), Combining Eqs. 3 , 4 , and 5 with Eq. 2 results in Setting V _{0} = 0 provides a maximum estimate of divergence time. Using V _{0} = (V _{1} + V _{2})/2 reduces Eq. 6 to a mutationdriftequilibrium model (Eq. 1 ), which provides a minimum estimate unless there was a decline in population size (23).
The rationale behind the three additional variance terms can be illustrated with an example. Imagine two daughter populations with very small effective population sizes that are at mutation drift equilibrium. Because of their small sizes, these populations will be very homogeneous. If they expand rapidly to large effective sizes, at first they will remain homogeneous, but over time, mutations will continue to introduce variation until a new equilibrium is reached. In the initial generations after the split, the morphological divergence will be very slow, because genetic drift will be very weak. This is because in a large population, the chance effects of sampling are small, and if the population is homogeneous, it does not matter much which individuals give rise to the next generation. As the populations become less homogeneous with the addition of mutational variance over time, the rate of divergence will increase until the mutationdriftequilibrium rate is reached. Because the divergence is slow initially, assuming that the populations were always at mutation, drift equilibrium will underestimate the actual divergence time. The three additional variance terms correct for this bias.
Effective Population Size.
If the amount of additive genetic variance introduced by mutation per generation, V _{m}, is known, then it should be possible to estimate the effective population size of a particular group (populations or species) from the amount of additive genetic variation found within that group. Let N _{e} be the effective size of a population. Following Lynch and Hill (20), at mutation drift equilibrium, the expected amount of additive genetic variance, V, found within the population is Eq. 7 has been tested empirically, and it appears to hold reasonably well, at least on average, in inbreeding experiments on Drosophila melanogaster (30). Substituting V for V _{1} and σ for σ_{1} in Eq. 3 , combining Eqs. 3 and 5 with Eq. 7 , and solving for N results in an estimator for the effective size of the population given by
Statistical Analyses
Data.
Our analyses are based on 37 standard cranial measurements (Fig. 1) collected on 2,524 modern humans from 30 globally distributed populations and 20 Neandertal specimens. More details about the sample and measurements can be found elsewhere (12, 31–33). We have previously shown that these 37 measurements appear to be evolving neutrally in Neandertals and modern humans (12), making them candidates for use in estimating divergence times.
Treatment of Multiple Measurements.
Divergence time can be estimated from multiple measurements as P̅T̅_{D}, the mean of PT _{D} estimates for individual measurements. The expected value for PT _{D} for each measurement is the divergence time, so the expected value for P̅T̅_{D} is also the divergence time, regardless of whether the measurements are correlated with each other. If the measurements are completely independent (uncorrelated within and between groups), bootstrapping (34) can be used to estimate the variance of P̅T̅_{D}, but this procedure will underestimate Var {P̅T̅_{D}} if the measurements are not independent.
It is always possible, however, to find a basis for which the measurements are uncorrelated within groups as the eigenvectors of the pooled within group variance–covariance matrix. There are as many eigenvectors as original measurements. Each eigenvector (sometimes called a principal component) represents a new measurement that is a linear combination of the original measurements and is uncorrelated with the other eigenvectors within groups. For neutral divergence, the elements of the betweengroup variance–covariance matrix are expected to be proportional to the elements of the withingroup variance–covariance matrix (35–37), so the eigenvectors should be approximately uncorrelated between groups as well. Eq. 8 for effective population size can be extended to multiple measurements in a similar manner.
Calibration.
As with divergence time estimates from molecular data, to estimate an “unknown” divergence time, we must first calibrate certain parameters using “known” reference points. In our case, we need to calibrate m and h ^{2}. We pick two “known” reference points. The first is P̅T̅_{D} under a mutationdriftequilibrium model for the split between subSaharan African and other human populations. This divergence is the oldest among human populations, so it is the most appropriate calibration point for estimating even older divergence times. Zhivotovsky (23) estimated T̅_{D} as ≈30,000 years ago based on 131 autosomal and four Ychromosomal microsatellites and from comparisons with other studies. This divergence time estimate is almost certainly too recent, because the mutationdriftequilibrium model assumes that there was no postdivergence gene flow between subSaharan African and other human populations, and that human populations have not expanded in size recently. However, by calibrating to a divergence time based on a mutationdriftequilibrium model, we are just assuming that the morphological and microsatellite estimates should match up under the same model, not that this is the most realistic model to use to infer the actual divergence time. In other words, because the mutationdriftequilibrium model is expected to underestimate the actual divergence time between subSaharan African and other human populations by the same amount for morphology and microsatellites, using 30,000 years ago as the calibration point will not result in an underestimate of the divergence time between Neandertals and modern humans.
The second reference point is the effective population size, P̅N̅_{e}, under a mutation–drift–equilibrium model for subSaharan African human populations. Zhivotovsky and colleagues (17) estimated N̅_{e} from 271 microsatellites using an equation equivalent to our Eq. 7 as ≈2,700 individuals. Once again, we are just assuming that the morphological and microsatellite estimates should match up under the same model, not that this is the most realistic model to use to infer the actual effective population size. Using V _{0} = (V _{1} + V _{2})/2 in Eq. 8 for mutation drift equilibrium and setting P̅T̅_{D} = T̅_{D} in Eq. 6 and P̅N̅_{e} = N̅_{e} in Eq. 8 produces two equations that can be solved for h ^{2} and m, the two remaining unknown parameters. Following Zhivotovsky (23), we use a generation length of 25 years for all our calculations.
Confidence Limits.
We use bootstrapping (34) to calculate approximate 95% confidence limits for P̅T̅_{D}. We resampled with replacement 10,000 times from the PT _{D} estimates for the individual eigenvectors and calculated an estimate of P̅T̅_{D} for each resample, producing a distribution of P̅T̅_{D} estimates. The lower and upper confidence limits are the 2.5% and 97.5% quantiles respectively of this distribution. These confidence limits account for evolutionary stochasticity in the amount of mutational variance actually introduced per generation, but they do not account for uncertainty in the calibration. Because uncertainty in the calibration is difficult to quantify accurately, it is not usually included in confidence limits for divergence times estimated from DNA sequences. We follow this practice here, but it is important to point out that this additional uncertainly can often be quite large.
Calculations.
We wrote scripts in MATLAB (Mathworks) to perform all of the calculations, making use of the statistical toolbox and Richard E. Strauss' MATLAB function package.
Results
Calibration.
From our calibration, we estimate h ^{2} = 0.37. Studies of human cranial measurements commonly use h ^{2} = 0.55 (5, 6, 38, 39), but this value actually derives from softtissue head measurements, which appear to have higher heritabilities than skeletal measurements (40). In a recent study of skeletal measurements collected on a pedigreed human sample, Carson (40) estimated heritabilities for 21 of the 37 measurements we use here with mean h ^{2} = 0.36. The sample sizes were small for three measurements, potentially making the h ^{2} estimates unreliable (40). For the remaining 18 measurements, mean h ^{2} = 0.31. Both estimates are quite close to our calibrationbased estimate, especially considering the considerable uncertainty in h ^{2} estimates. Additionally, we do not expect an exact match, because the populations are different, and our calibrationbased estimate reflects average h ^{2} over hundreds to thousands of generations (evolutionary time) rather than over just a few generations.
From our calibration, we estimate m = 1.20 × 10^{−4}. Experimental studies of a number of different taxa for a variety of traits (21, 29, 41) suggest a rough range of 10^{−4} to 10^{−2} for m. Our calibrationbased estimate is at the lower end of the range, which is expected, because most experimental values for m are overestimates. Experimental estimates consider all mutations, not just neutral ones, and over evolutionary time stabilizing selection will remove some mutations from the population. The overestimation may sometimes be substantial, because the experimental populations have low effective sizes, which would weaken stabilizing selection (1). A similar overestimation problem is documented for molecular evolution where mutation rates estimated from pedigrees are much higher than those estimated from phylogenetic calibration points (42, 43).
Neandertal and Modern Human Divergence.
We estimate that Neandertals and modern humans diverged ≈311,000 years ago (95% C.I.: 182,000–466,000) assuming mutation drift equilibrium or 435,000 years ago (95% C.I.: 308,000–592,000) assuming V _{0} = 0. For both estimates, we added 25,000 years to account for the fact (averaging dates) that Neandertals lived ≈50,000 years ago. When we compare Neandertals with only male recent humans, the point estimates and C.I.s decrease by <10,000 years, so our estimates would not be strongly biased, even if the entire Neandertal sample were male.
It is difficult to decide which V _{0} model is most appropriate. The V _{0} = 0 result is probably an overestimate for at least two reasons. (i) Because the Neandertal sample is too small to accurately estimate withinpopulation variation in Neandertals, we used the human value for both V _{1} and V _{2}. If Neandertals were actually less variable than presentday human populations, the V _{0} = 0 result would be an overestimate (i.e., the Neandertal lineage would maximally deviate from mutation drift equilibrium less than the modern human lineage). (ii) The additive genetic variance in the last common ancestor of Neandertals and modern humans must have been greater than zero, making the V _{0} = 0 result an overestimate. In contrast, the mutationdrift equilibrium, V _{0} = (V _{1} + V _{2})/2, result could be an underestimate for at least two reasons. (i) Human populations have grown in size recently, which would make the mutationdriftequilibrium result an underestimate as long as this growth in census size corresponds to growth in effective size. (ii) Postdivergence gene flow between Neandertals and modern humans would make the mutationdriftequilibrium result an underestimate (23). Given this uncertainty, the mean of the two estimates, 373,000 years ago, seems to be a reasonable point estimate.
Discussion
If we consider the maximum extent of the 95% confidence limits for both the mutationdrift equilibrium and the V _{0} = 0 estimates, then the Neandertal and modern human lineages split between 182,000 and 592,000 years ago. Although this range is quite large, it still allows for some observations with respect to the human fossil record. First, even the lower limit is within the Middle Pleistocene, suggesting a relatively deep divergence of Neandertals and modern humans, which is consistent with the presence of derived Neandertal features on Middle Pleistocene fossils from Europe (44–47). Second, recent dates suggesting that the Sima de los Huesos site is >530,000 years old (48) would put the fossils from this site, which appear to have multiple derived Neandertal features (46), at or potentially before the split of the Neandertal and modern human lineages. Third, although the ≈800,000yearold AtapuercaTD6 humans (49) could be ancestral to Neandertals and modern humans, their morphology may not be representative of the source population that actually gave rise to Neandertals and modern humans, because they date from, at minimum, ≈200,000 years before the split time.
Our PT _{D} results are estimates of when the ancestral Neandertal and modern human populations last shared a randomly mating common ancestor (split time), whereas most molecular estimates are DNA sequence coalescence times. As long as the mutation rate is correct, coalescence times will equal or predate the split time by an amount that depends, on average, on ancestral effective population size (50). For example, if the ancestral population had a constant effective size of 2,500 individuals, an autosomal coalescence time based on one Neandertal and one extant human sequence would be expected to predate the split time by 125,000 years, assuming a generation length of 25 years. Point estimates for the coalescence of ancient Neandertal and extant human sequences for both mitochondrial and nuclear DNA range from ≈300,000 to 800,000 years ago (42, 51–55). Additionally, Noonan and colleagues (53) estimated that the Neandertal and modern human lineages split ≈370,000 years ago based on comparisons at >36,000 autosomal DNA sites. This split time is very similar to 373,000 years ago, the average of our mutationdriftequilibrium and V _{0} = 0 point estimates.
The divergence time estimates could change somewhat if the separation between Neandertals and modern humans was actually more complicated than a simple splitting of populations. Additionally, although the ancestors of Neandertals and modern humans may have split ≈400,000 years ago, modern humans are sometimes thought to have originated with a subsequent speciation event ≈150,000 years ago. It is unclear what this view implies demographically, but it could be taken to mean that the emergence of modern humans involved a bottleneck, which is a sharp reduction followed by a rapid expansion in population size. If such a bottleneck occurred, the mutationdriftequilibrium result would likely be an underestimate, and the V _{0} = 0 estimate would be closer to the actual split time.
Regardless of these potential complications, the close correspondence between the cranial and DNAsequence estimates implies that both datasets largely, although not necessarily exclusively, reflect neutral divergence, causing them to track population history or phylogeny rather than the action of diversifying natural selection. The overall pattern for the measurements considered here appears to be neutral divergence, but as is the case among human populations (6–8, 38), a few measurements could still have diverged by diversifying natural selection. Additionally, future research will be needed to determine if this pattern is representative of Neandertal and modern human cranial divergence in general, or whether it applies only to the measurements included in our study. In particular, more detailed or internal cranial structures and other aspects of cranial anatomy that are difficult to quantify with standard osteometric tools many have been shaped by diversifying natural selection.
Brain size relative to body size appears to have increased in both the Neandertal and the modern human lineages (56, 57). These parallel trajectories may indicate that directional natural selection was acting on both lineages independently, resulting in differently shaped but similarly sized brain cases (58). To the extent that the measurements considered here reflect these changes, our results would imply that although natural selection may have produced the similarities in size, genetic drift lead to the differences in shape. We are not arguing that natural selection had no effect, just that it appears not to have played a dominant role in producing differences.
One misconception about neutral evolution is that it is slow, because it is driven by genetic drift rather than by natural selection. However, if stabilizing natural selection is more prevalent than directional natural selection, then neutral evolution will actually be comparatively fast. Lynch (59) estimated that human cranial evolution was rapid relative to morphological evolution in other mammals, but the absolute rate was consistent with neutral evolution. It follows from these results and ours that Neandertal and modern human crania may have been released from the constraints of stabilizing natural selection that limit the rates of morphological evolution in other mammals.
Perhaps the most significant implication of our results is that there is no conflict between molecules and morphology. This contrasts with a common characterization of debates about the origins of modern humans as molecules vs. morphology (60). In fact, at least for the measurements considered here, there is a close quantitative correspondence between the amount of cranial divergence and the amount of DNA sequence divergence between Neandertals and modern humans.
Acknowledgments
We thank D. DeGusta, R. Fournier, B. Henn, J.J. Hublin, K. A. Horsburg, R. Klein, J. Mountain, N. Rosenberg, T. Steele, M. Stoneking, and two anonymous reviewers for valuable comments on drafts or feedback; the late W. W. Howells for making his raw data publicly available; the institutes that made their fossil material available for study; M. v. Harling for the CT scans used to create Fig. 1; and the Max Planck Society for financial support.
Footnotes
 ^{‡}To whom correspondence should be addressed. Email: tdweaver{at}ucdavis.edu

Author contributions: T.D.W., C.C.R., and C.B.S. designed research; T.D.W. and C.C.R. performed research; C.B.S. contributed new reagents/analytic tools; and T.D.W., C.C.R., and C.B.S. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

↵ ‖ Because of a spelling change in German, there are two acceptable spellings for the name of the group. “Neanderthal” rather than “Neandertal” is the preferred spelling of C.B.S.
 © 2008 by The National Academy of Sciences of the USA
References
 ↵
 ↵
 ↵

↵
 Lynch M
 ↵

↵
 Roseman CC
 ↵

↵
 Harvati K ,
 Weaver TD
 Harvati K ,
 Harrison T
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Ramachandran S ,
 et al.
 ↵

↵
 Lande R

↵
 Lande R

↵
 Lynch M ,
 Hill WG

↵
 Turelli M ,
 Gillespie JH ,
 Lande R

↵
 Goldstein DB ,
 Ruiz Linares A ,
 CavalliSforza LL ,
 Feldman MW

↵
 Zhivotovsky LA

↵
 Zhivotovsky LA ,
 Feldman MW

↵
 Goldstein DB ,
 Linares AR ,
 CavalliSforza LL ,
 Feldman MW

↵
 Slatkin M

↵
 Kremer A ,
 Zanetto A ,
 Ducousso A

↵
 Falconer DS
 ↵

↵
 Whitlock MC ,
 Fowler K

↵
 Howells WW

↵
 Howells WW

↵
 Howells WW
 ↵

↵
 Lande R
 ↵

↵
 Ackermann RR ,
 Cheverud JM
 ↵
 ↵
 ↵

↵
 Houle D ,
 Morikawa B ,
 Lynch M

↵
 Ho SYW ,
 Phillips MJ ,
 Cooper A ,
 Drummond AJ
 ↵

↵
 Hublin JJ
 Akazawa T ,
 Aoki K ,
 BarYosef O
 ↵
 ↵
 ↵
 ↵

↵
 Bermúdez de Castro JM ,
 et al.

↵
 Rosenberg NA ,
 Feldman MW
 Slatkin M ,
 Veuille M
 ↵
 ↵

↵
 Noonan JP ,
 et al.
 ↵

↵
 Beerli P ,
 Edwards SV
 ↵
 ↵

↵
 Bruner E ,
 Manzi G ,
 Arsuaga JL
 ↵
 ↵
Citation Manager Formats
Article Classifications
 Biological Sciences
 Anthropology