Placing confidence limits on the molecular age of the human–chimpanzee divergence

Results and Discussion

We first used the ML-distance method in which a lineage-specific molecular clock is applied to estimate the human– chimpanzee divergence time (see Fig. 1 and Materials and Methods). This method requires a priori knowledge of at least one species divergence time to calibrate the local molecular clock. We took a conservative approach and used a minimum date of 23.8 Ma (boundary of the Oligocene and Miocene) for the ape–OWM divergence (28) (see Supporting Text). We focused on the minimum time of divergence because the inference of the true (mean) time of divergence of any two lineages requires a more robust fossil record for calibration than is generally available for most groups of organisms, and especially these primates. Most past molecular clock studies have also estimated minimum divergence times for this reason (29). In recent years, many methods have been developed that attempt to estimate mean times of divergence by specifying both minimum and maximum calibration points. However, whereas most minimum calibrations are robust, this is not true for maximum calibrations or the implicit assumption about the probability density assumed for fossil calibration uncertainty (29). It is difficult to provide convincing evidence that a phylogenetic divergence event did not occur earlier than a particular point in time. Errors in assignment of maximum constraints, whether too young or too old, will potentially bias the mean time estimate.

Analysis of Third Codon Positions. A likelihood ratio analysis of third codon positions from 167 genes (53,008 codons) indicated significant rate variation among sites (Fig. 1e ), so we used a GTR+Γ model in all further analyses. Using a 23.8 Ma minimum date for ape–OWM divergence, we obtained an ML-distance estimate of 4.74 Ma, which is older than the age of 4.2 Ma for the earliest unequivocal postcranial evidence of upright hominids (14) but younger than some recently proposed early hominids dated from between 5 and ≈6 Ma (15–19).

However, the absolute dates derived above depend directly on the time used to calibrate the molecular clock, and it is important to focus on the size of the C.I. relative to the estimate, rather than on the point estimate of time alone. Using the MBR approach, we obtain a 95% C.I. of 3.88–5.65 Ma for the minimum date. This C.I. spans –18% to +19% of the estimated time. Our computer simulations under equal, random, and correlated models of rate evolution show that the 95% C.I. generated by MBR indeed contains the true value ≥95% of the time when all third positions are resampled, regardless of gene boundaries (Fig. 2).

Of all third codon positions, the 4-fold-degenerate sites are expected to be under the least amount of natural selection. We divided these sites into those that were involved in CpG dinucleotides and those that were not, because the 4-fold-degenerate sites involved in CpG dinucleotide configurations mutate much faster (24, 25). Analysis of 19,311 non-CpG 4-fold-degenerate sites produced an estimate of 4.75 Ma, which is almost identical to that given by third codon positions, but the 95% C.I. was wider (–25% to +31%). Because the data set size (in terms of the number of sites) was the largest for the third codon positions, we considered only third codon positions in all further DNA sequence analyses.

Analysis of Nonsynonymous Substitutions. Unexpectedly, ML-distance analyses of amino acids and second codon positions yielded estimates of human–chimpanzee divergence that were >40% larger than those obtained from the third codon positions and 4-fold-degenerate sites. In fact, neither the estimate obtained using amino acids (6.80 Ma), nor the third-position estimate (4.74 Ma) fell within the 95% confidence limit of the other (4.78–8.86 Ma and 3.88–5.65 Ma for amino acids and third positions, respectively). This finding is surprising because both DNA and protein substitutions are expected to provide concordant results, especially when using a large number of protein-coding genes.

A cursory examination of the time estimates based on the rate of amino acid substitutions revealed a tendency toward larger estimates for very slowly evolving proteins. The 40 most slowly evolving proteins (based on evolutionary divergence of human and mouse proteins) provide an estimate more than three times higher than the 40 fastest evolving proteins. In contrast, the third codon position estimates are similar whether obtained from genes corresponding to proteins evolving with low, medium, or fast rates and coincide with the estimates obtained from fast-evolving proteins (Fig. 3). These patterns suggest a disproportionate inflation of sequence divergence for closely related species as compared with distantly related species, an effect that would be the greatest for the most slowly evolving proteins. The skewed inflation of sequence divergence estimates for closely related sequences may be attributed to sequencing errors, presence of ancestral polymorphisms (30), persistence of slightly deleterious polymorphisms (31, 32), and masking of deleterious mutations in heterozygotes (33). Although these factors affect all amino acids and third codon positions, the most slowly evolving sequences (e.g., highly conserved proteins) are affected the most relative to their rate and the time of sequence divergence.

Our computer simulations showed that if the cumulative effect of nonneutral polymorphism due to factors described above produced, at an average, three differences per 1,000 aa, then an overestimation pattern similar to the magnitude and trend observed for the real data would be obtained (Fig. 3, gray area). Furey et al. (34) have already reported sequencing error rates of 0.5 bp per 1,000 nucleotides for human GenBank data, which translates to one error per 1,000 aa (because >75% of base pair errors in three codon positions result in an error at the amino acid sequence level). Although genomic sequences have a lower error rate, within-species comparisons in humans indicate a 1 in 1,000 bp difference at the genome level between individuals (35). This variation contributes nonfixed differences between species and affects the estimate of human–chimpanzee divergence in the same way as the sequencing errors. Sequencing error rate for available chimpanzee sequences is higher than that for humans (4 times versus 10 times coverage; see Ensembl, www.ensembl.org), and chimpanzees are known to have a much higher within-population polymorphism (33), both of which would lead to an overestimation in the present case.

In addition to the above factors, the fixation of lineage-specific amino acid substitutions due to positive selection (36) may also make the use of amino acid sequences problematic. However, it is unlikely to be the major contributing factor in the present case, because none of the slowest-evolving proteins in our data set were found to be involved in host defense, immunity, or olfaction (36, 37). On the other hand, the most rapidly evolving proteins, which show estimates similar to those obtained from the third codon positions, contain a large proportion of genes known to be candidates for positive selection (which show nonsynonymous to synonymous substitution ratios of >0.5). In any event, because of the observed artifact mentioned above, we do not further consider the estimates based on amino acid sequence analysis.

Bayesian Inference of the Minimum Time for Human–Chimpanzee Divergence. The relative rate analysis shows that human third codon positions are 17% more divergent than those of chimpanzee and that the 4-fold-degenerate non-CpG sites are evolving ≈20% more slowly in hominoids than in OWM. Both of these rate differences are insignificant [P > 0.05 in Tajima's test (23)] and are in agreement with those reported elsewhere (33). The difference between hominoids and OWM is similar in magnitude to that reported in an analysis of large data sets by Yi et al. (38) and Kumar and Subramanian (39). It is one-half of that reported by Steiper et al. (40), whose estimates of relative rates depend on fossil-based divergence times for within-hominoid and within-cercopithecoid species, although an outgroup was not used in their analyses.

In any case, the observation of any rate differences in the hominoid–OWM comparison means that the ancestral human–chimpanzee lineage may not be evolving at the same rate as those leading to humans and chimpanzees, violating the assumptions made in the ML-distance approach. Therefore, we applied Bayesian analysis as implemented in multidivtime (21) that models rate variation among lineages to generate a better estimate of human–chimpanzee divergence. We used 23.8 Ma for the ape–OWM divergence time as the root-to-tip mean (RTTM = 23.8) in the multidivtime Bayesian analyses to make results directly comparable to the ML-distance method (see Supporting Text for parameters used for Bayesian analysis). This produced an estimate of 4.98 Ma for the third codon positions and of 5.17 for the 4-fold-degenerate non-CpG sites. These estimates are marginally higher than those for the ML-distance estimates, but the discrepancy becomes smaller when we compare the ratio of the times of the target (human–chimpanzee) and calibration (ape–OWM) divergence events. Even though we provide an explicit time (RTTM) for the ape–OWM split, multidivtime infers time estimates for both human–chimpanzee and ape–OWM divergences. These estimates give consistent ratios of 0.21 (4.98/24.00) and 0.21 (5.17/24.37) for the third codon positions and the 4-fold-degenerate non-CpG sites, respectively. In the ML-distance method, the time ratio is 0.20 for the third codon positions (4.74/23.8) and also 0.20 (4.75/23.8) for the 4-fold-degenerate non-CpG sites. Therefore, the time ratios for the Bayesian estimates are only 5% higher than are those obtained using the ML approach. Bayesian methods also produced a narrower 95% C.I. (–12% to +19%) as compared with the C.I. (–18% to +19%) of the ML-distance approach for the third codon position data. Our computer simulations confirmed that the MBR C.I.s contain the true value with a frequency of >95% in the Bayesian approach as well (Fig. 2). Because the C.I.s generated using the MBR approach with the Bayesian analysis are smaller, we use them below.

Comparison with Previous Studies. Our minimum estimates of 4.74 and 4.98 Ma from ML-distance and Bayesian methods, respectively, based on the minimum 23.8 Ma date for ape–OWM divergence, are smaller than many previous molecular time estimates. There are three primary reasons for this difference. First, many studies have presented (and sometimes preferred) the higher divergence times estimated by using the protein sequences or by using all three codon positions (e.g., refs. 5 and 10). However, as shown above, the inclusion of slowly evolving amino acid sequences (or first and second codon positions) is likely to bias time estimates for the human–chimpanzee divergence. This pattern is also seen in the Bayesian analysis of amino acid sequences, which produces an ≈60% larger time estimate as compared with the third codon positions, with a 95% C.I. that does not overlap with that obtained from the analysis of third codon positions.

Second, different studies have used vastly different divergence times (20–35 Ma) for the ape–OWM split to calibrate molecular clocks. A comparison of the estimated and assumed time ratios resolves this discrepancy. For example, our ratios from 0.21 to 0.20 are consistent with the 0.21 value reported in analyses conducted by using >150,000 bp of noncoding data by Stieper et al. (40), 97 protein-coding genes by Wildman et al. (9), and complete mitochondrial DNA by Schrago et al. (ref. 41; see also ref. 7). Therefore, in the absence of a firmly established calibration date, it is better to consider all results as ratios of human–chimpanzee to ape–OWM divergence times.

Third, most previous studies were based on the analysis of fewer genes. In those cases, the true C.I. accounting for the variance from gene sampling, multiple-hit correction, and evolutionary rate differences among species is expected to be rather large. This effect is evident from the lower and upper limits of the 95% C.I.s obtained by using smaller random subsets of 167 genes (Fig. 4). When the number of genes is <10 (e.g., ref. 5), the 95% C.I. is nearly as large as the time estimate itself. The C.I. declines to ≈50% if 50 genes are used (e.g., refs. 7 and 10). With the use of 100 or more genes, uncertainties contributed by rate variation among lineages become predominant, and the C.I. width declines only slowly. Thus, little is to be gained by adding data from more genes, but using more species may help narrow the 95% C.I. width in the future.

Incorporating Multiple Fossil Bounds in the Bayesian Inference. Bayesian methods allow more fossil constraints to be used in species divergence time estimation, such as a minimum constraint of 4.2 Ma (14) for the human–chimpanzee divergence. After adding this constraint and keeping RTTM = 23.8, we found that the human–chimpanzee divergence time increased by 53%. However, the ratio of estimated times for ape–OWM and human–chimpanzee divergences was 0.21, which is the same ratio obtained without placing the minimum bound on the human–chimpanzee divergence. Setting RTTM (time for the common ancestor of human, chimpanzee, and macaque) to 35 Ma in the Bayesian analyses again produced the same time ratio (0.21), even though the estimate of human–chimpanzee time became even larger.

To determine which estimates are likely to be correct, we conducted analyses with and without the human–chimpanzee constraint of 4.2 Ma for the simulated data under equal, random, and correlated evolutionary rate variation. [In the random rate case, the rate at each branch is randomly selected from a uniform distribution of rates so as to introduce a ±40% random noise in evolutionary rate independently for each gene. In the correlated rate, the assignment of lineage-specific rate uses a stationary lognormal distribution in which the rates vary from branch to branch as a random walk for a given gene, so that rates drift up or down along the branches of any lineage (42, 43). See Supporting Text.] In all cases, time ratios were recovered with high accuracy, but the absolute times estimated were strongly dependent on the calibration constraints used in the Bayesian analyses, as observed in the analysis of the real data.

At this stage, it is important to revisit the distinction between estimating minimum and true times of species divergence (see discussions in ref. 29). Clearly, the exclusive use of minimum calibration times, derived from the fossil record, will produce only minimum time estimates of the target divergence and its C.I., irrespective of the sophistication of the methods used (29). Focusing on minimum time estimates provides at least two advantages. For one, asserting a minimum time constitutes a definite, falsifiable statement, namely that a divergence occurred before some specified date. Alternatively, we could specify upper and lower bounds on calibration times (44) in an effort to estimate “true” divergence time, but doing so would require the incorporation of invariably subjective specifications of calibration time uncertainty. The results of our application of Bayesian methodologies to real and simulated data, discussed above, indicate that the absolute times of human–chimpanzee divergence inferred by using a four-species construct are strongly affected by the upper and lower bounds used, but these methods are successful in correctly estimating ratios of time estimates for different nodes in the phylogeny.

In fact, the estimated C.I.s from Bayesian analyses are unlikely to be around the “true” divergence time, unless a probability distribution around the range of times inferred from the fossil record can be assigned (29). For example, if we have established minimum and maximum values for calibration time but have reason to believe from fossil evidence that the true divergence time is closer to the former than the latter, then the time distribution may be triangular or lognormal (29). In other cases, the best we can do is to associate a standard deviation with the calibration time and, if necessary, assume a normal or uniform distribution. In some instances we may have evidence pointing toward a certain shape of distribution, in which case the MBR method can be used to explicitly specify this distribution. For example, if we assume a linearly declining probability distribution with 23.8 Ma as minimum (and most probable) estimate of ape–OWM and 35 Ma as its maximum, we obtain a 95% C.I. of 4.17–7.13 Ma. In this case, the C.I., rather than the actual date, should be reported because the ad hoc form of distribution used strongly determines the central value of the time estimated rather than any biological knowledge.

Conclusions

Even when mean or true times of divergence between species are difficult to estimate with molecular clocks, because of uncertainties in fossil calibrations, minimum and relative times of divergence can still be obtained. Here we have shown that the divergence of chimpanzees and humans is such a case, because of the current lack of robust calibration points. Our genome-scale analyses show that the minimum time of that divergence was very close to one-fifth of the ape–OWM divergence, with a 95% C.I. from –12% to +19% of that time. This result can be used to test hypotheses that arise in the literature. For example, if claims of late Miocene (≈6 Ma) fossil hominids (15–17) are correct, then the ape–OWM divergence occurred at least 24–25 Ma or earlier based on the 95% C.I. (Table 1). This date extends the age of that latter divergence only slightly from the accepted 23.8 Ma. Conversely, if fossils of hominoids or OWM were discovered and dated to 27 Ma, this would imply a minimum time of divergence of between 4.9 and 6.6 Ma for chimpanzees and humans. Because large gaps in the fossil record are not implied in either case, they would be compatible with current understanding of the primate fossil record. In contrast, the recent suggestion that some ape-like fossil teeth from Africa might be from the African ape clade and that the human–chimpanzee divergence time might be even older than 12.5 Ma (20) seems to be extremely unlikely. This inference is because extrapolation of the upper 95% confidence limit from Table 1, using this date, implies a minimum ape–OWM divergence time of ≈55 Ma, >30 million years earlier than supported by current fossil evidence. Similarly old molecular clock dates (13.5 Ma, based on mitochondrial DNA and using nonprimate calibration points) for the human–chimpanzee divergence (12, 45) present the same problem. By enabling tests such as these, new fossil discoveries and better radiometric dating of existing fossils provide the best prospects for improved understanding of the timing of events in hominid evolution.