Edited by Erik Trinkaus, Washington University, St. Louis, MO, and approved January 17, 2008 (received for review September 27, 2007)
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 within-population variation. These dates are quite similar to those recently derived from ancient Neandertal and extant human DNA sequences. Close correspondence between cranial and DNA-sequence 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 among-population 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, within-population 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.
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 between-population 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 between-population 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 within-population 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 among-individual 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 narrow-sense 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 mutation-drift-equilibrium 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 mutation-drift-equilibrium 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.
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
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.
The approximate locations of the cranial measurements used in the analyses are superimposed as red lines on lateral (A), anterior (B), and inferior (C) views of a human cranium. Note that when the endpoints of a measurement are not visible, the line is projected into a plane situated in front of the cranium.
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 between-group variance–covariance matrix are expected to be proportional to the elements of the within-group 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.
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 mutation-drift-equilibrium 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 Y-chromosomal microsatellites and from comparisons with other studies. This divergence time estimate is almost certainly too recent, because the mutation-drift-equilibrium 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 mutation-drift-equilibrium 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 mutation-drift-equilibrium 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.
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.
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.
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 soft-tissue 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 calibration-based 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 calibration-based 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 calibration-based 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).
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 within-population variation in Neandertals, we used the human value for both V 1 and V 2. If Neandertals were actually less variable than present-day 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 mutation-drift 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 mutation-drift-equilibrium 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 mutation-drift-equilibrium result an underestimate (23). Given this uncertainty, the mean of the two estimates, 373,000 years ago, seems to be a reasonable point estimate.
If we consider the maximum extent of the 95% confidence limits for both the mutation-drift 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,000-year-old Atapuerca-TD6 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 mutation-drift-equilibrium 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 mutation-drift-equilibrium 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 DNA-sequence 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.
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.
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.
