Strong mitochondrial DNA support for a Cretaceous origin of modern avian lineages
© Brown et al. 2008
Received: 02 January 2008
Accepted: 28 January 2008
Published: 28 January 2008
Skip to main content
© Brown et al. 2008
Received: 02 January 2008
Accepted: 28 January 2008
Published: 28 January 2008
Determining an absolute timescale for avian evolutionary history has proven contentious. The two sources of information available, paleontological data and inference from extant molecular genetic sequences (colloquially, 'rocks' and 'clocks'), have appeared irreconcilable; the fossil record supports a Cenozoic origin for most modern lineages, whereas molecular genetic estimates suggest that these same lineages originated deep within the Cretaceous and survived the K-Pg (Cretaceous-Paleogene; formerly Cretaceous-Tertiary or K-T) mass-extinction event. These two sources of data therefore appear to support fundamentally different models of avian evolution. The paradox has been speculated to reflect deficiencies in the fossil record, unrecognized biases in the treatment of genetic data or both. Here we attempt to explore uncertainty and limit bias entering into molecular divergence time estimates through: (i) improved taxon (n = 135) and character (n = 4594 bp mtDNA) sampling; (ii) inclusion of multiple cladistically tested internal fossil calibration points (n = 18); (iii) correction for lineage-specific rate heterogeneity using a variety of methods (n = 5); (iv) accommodation of uncertainty in tree topology; and (v) testing for possible effects of episodic evolution.
The various 'relaxed clock' methods all indicate that the major (basal) lineages of modern birds originated deep within the Cretaceous, although temporal intraordinal diversification patterns differ across methods. We find that topological uncertainty had a systematic but minor influence on date estimates for the origins of major clades, and Bayesian analyses assuming fixed topologies deliver similar results to analyses with unconstrained topologies. We also find that, contrary to expectation, rates of substitution are not autocorrelated across the tree in an ancestor-descendent fashion. Finally, we find no signature of episodic molecular evolution related to either speciation events or the K-Pg boundary that could systematically mislead inferences from genetic data.
The 'rock-clock' gap has been interpreted by some to be a result of the vagaries of molecular genetic divergence time estimates. However, despite measures to explore different forms of uncertainty in several key parameters, we fail to reconcile molecular genetic divergence time estimates with dates taken from the fossil record; instead, we find strong support for an ancient origin of modern bird lineages, with many extant orders and families arising in the mid-Cretaceous, consistent with previous molecular estimates. Although there is ample room for improvement on both sides of the 'rock-clock' divide (e.g. accounting for 'ghost' lineages in the fossil record and developing more realistic models of rate evolution for molecular genetic sequences), the consistent and conspicuous disagreement between these two sources of data more likely reflects a genuine difference between estimated ages of (i) stem-group origins and (ii) crown-group morphological diversifications, respectively. Further progress on this problem will benefit from greater communication between paleontologists and molecular phylogeneticists in accounting for error in avian lineage age estimates.
Many evolutionary models [1–4] are tightly linked to absolute timescales. A reliable temporal framework is therefore required for understanding the tempo (and, through correlation with geophysical phenomena, mechanisms) of biological evolution. There are two complementary sources of information for dating ancient biological divergences: (1) physical historical remains (either paleontological or ichnological); and (2) molecular sequence differences among extant taxa, the analysis of which requires assumptions about the processes and rates of sequence evolution. Unfortunately, these two sources of information ('rocks' and 'clocks', respectively) often yield starkly disparate estimates of the timing of major biological divergences .
Strictly speaking, any molecular estimate that generates a positive value of δ Realized MA-FA is consistent with the fossil record. It is instead the magnitude of δ Realized MA-FA that suggests conflict, and distressingly large values sometimes exist. Conflict between the two sources of information is especially evident with respect to the ages of extant avian lineages (Neornithes). Based on a strict interpretation of the fossil record (i.e. insignificant δ Fossil error ), Feduccia [8, 9] proposed an explosive Cenozoic origin for most modern avian lineages, presumably a result of open niches left by newly extinct non-avian dinosaurs and other taxa. Although a recent fossil find  forces a minimum of five of the earliest Neornithes divergences into the late Cretaceous, the fossil record generally supports the view that most modern lineages originated in the Cenozoic [8, 9, 11–15]. In contrast, molecular estimates all indicate that these same lineages are considerably older, sometimes as much as twice as old as analogous paleontological estimates [4, 16–26]. Between these two extremes lies the Cretaceous-Paleogene (K-Pg; formerly Cretaceous-Tertiary or K-T) boundary, a period when as many as 50% of land-dwelling species went extinct . The conflicting age estimates thus have different implications regarding the influence of the K-Pg mass extinctions on the evolutionary radiation of modern birds.
Although resolution of this conflict is clearly important for understanding avian diversification, it is hindered by compelling arguments from both sides. The supposition that the quality of the fossil record deteriorates backwards in time, for example, is contradicted by the observed congruence between stratigraphic and phylogenetic ordering of taxa . Sophisticated stratigraphic analyses indicate that fossils of the antiquity necessary to produce congruence with molecular studies are extremely improbable [11, 29, 30] (but see [31, 32]). Furthermore, of the known Mesozoic avian fossils [12, 14, 33, 34], the vast majority unambiguously lay outside Neornithes . Although a few Cretaceous fossils putatively represent modern lineages (e.g. parrot , loon  and others [12, 14]) these have largely been dismissed as fragmentary and inconclusive [9, 12, 38, 39]. One the molecular side, the recognition that rates of molecular evolution are often not clock-like (including birds [23, 40–42]), and that lineage-specific heterogeneity is common , has spurred the development of numerous 'relaxed' molecular clock methods (see reviews in [44–46]). In support of molecular genetic data, these methods perform well in simulation [47, 48] and, when applied to empirical data, deliver Cretaceous ages for the origin of modern birds [16, 23].
Given these arguments, the paleontological and molecular phylogenetic communities are currently at an impasse regarding the application of an absolute temporal axis for early organismal evolution [33, 49], and a range of evolutionary models [1–4] remain viable for birds. Here we endeavour to determine whether a more rigorous treatment of molecular genetic data lessens the 'rock-clock' discrepancy (δ Realized MA-FA ). In particular, we incorporate large fossil and taxon data sets, two components of molecular dating that have been shown to have a strong impact on the resulting divergence time estimates [50, 51]. In addition, we accommodate and explore the impact of uncertainty in both tree topology and molecular dating strategy. Finally, we test for signals of episodic molecular evolution related to both speciation events and absolute geologic time, processes that could potentially mislead molecular-based age estimates by systematically inflating branch lengths within speciose clades .
Degree of autocorrelation in rates of molecular evolution by partition and tree topology as calculated in Multidivtime
Autocorrelation (95% CI)
0.00197 (0.00127, 0.00290)
0.00437 (0.00258, 0.00685)
0.00452 (0.00288, 0.00680)
0.00566 (0.00343, 0.00874)
0.00177 (0.00112, 0.00263)
0.00344 (0.00197, 0.00548)
0.00380 (0.00241, 0.00571)
0.00414 (0.00206, 0.00744)
Model comparisons for analyses relaxing the assumption of autocorrelation of rates across the tree. Harmonic mean model likelihoods were calculated from post-burnin MCMC samples generated in BEAST. For these model comparisons, the topology was fixed as T Consensus. The strict clock model serves as a base comparison. The tree T Flexible refers to analyses where topology is not fixed. Covariance indicates the degree of substitution rate autocorrelation between ancestor and descendent branches. 95% HPDs are given in parentheses.
0.039 (-0.103, 0.175)
0.066 (-0.061, 0.193)
0.042 (-0.071, 0.161)
Estimated divergence times for major avian clades compared across methods and topologies. Estimated time to most recent common ancestor (tMRCA) of clades of non-controversial monophyletic status. Standard deviations are given in parentheses (for Dating5 and BEAST, standard deviations were calculated from 95% confidence/credibility intervals using a normal approximation). For r8s, PATHd8 and Multidivtime ages were estimated on each of the two fixed topologies (T Consensus and T Optimal). For BEAST, divergence times were estimated simultaneously with phylogeny (T Flexible). For each method an estimate of the magnitude of disagreement between fossil and molecular estimates of divergence times (δ Realized MA-FA ) is calculated as an average of MA-FA (the molecular age minus the fossil age) for those 18 internal nodes with calibration points.
δ Realized MA-FA
Second, the choice of the relaxed clock method had a strong influence on inferred ages. R8s, Multidivtime and BEAST tended to deliver similar estimates for most clades of interest (Table 3). In contrast, PATHd8 generated considerably younger dates with much smaller confidence intervals, despite using the same bootstrapped phylograms and fossil constraints as r8s. Dating5 tended to produce the most extreme results, with inferred basal split estimates similar to those from Multidivtime, but some derived split estimates younger than those from PATHd8. Most significantly, PATHd8 and Dating5 together identified five of these major clades as having crown diversification restricted to the Cenozoic (Ratites, Charadriiformes, Procellariiformes, Cuculiformes and Apodiformes), although the remaining methods generate estimates for these same nodes that are on average more than 50% older. In terms of comparing molecular and fossil age estimates (average δ Realized MA-FA ), r8s, Multidivtime and BEAST all show considerable discordance between the two sources of data, with the average molecular estimates for the major nodes (Table 3) being 36–45 MY older than corresponding fossil ages. PATHd8 and Dating5, in contrast, exhibit greater agreement between estimates from 'rocks' and 'clocks', with an average discrepancy of 17–25 MY.
Whether using fossil or molecular data, a daunting impediment to divergence time estimation in birds is that resolution of many inter-ordinal phylogenetic relationships has proven obstinate, despite large data matrices with multiple character types . Although our reconstruction T Optimal is overly optimistic in being fully resolved, it provides a useful alternative to the conservative T Consensus (Figure 2).
TOptimal recovers several traditional orders as polyphyletic (Caprimulgiformes, Coraciiformes, Falconiformes, Ciconiiformes), consistent with expectations  (but see ). Although TOptimal has caprimulgiform (nightbirds) families much more distantly related to one another than previous morphological  and molecular  investigations, differences in taxon sampling confounds direct comparison across studies. While Coraciiformes (kingfishers and relatives) is not found to be monophyletic, the two recovered sub-groupings both fall within a larger clade containing owls (Strigiformes), parrots (Psittaciformes) and woodpeckers and relatives (Piciformes). The monophyletic status of the order Falconiformes has received mixed support in previous analyses [54, 57–61]. Placement of Falconidae in TOptimal is suspect and likely stems from insufficient taxon sampling from this family . Regardless, no support was found for hypotheses uniting falconiform taxa with owls (Strigiformes)  or New World vultures (Cathartidae) with storks (Ciconiiformes) .
Several monotypic avian families have traditionally proved difficult to classify. The enigmatic hoatzin (Opisthocomidae) has variously been allied with at least four distantly related orders (Galliformes, Cuculiformes, Musophagiformes and Columbiformes; see [61, 63]). We find here an alliance with doves (Columbiformes), similar to joint analyses of mitochondrial and nuclear DNA sequences . The taxonomically problematic sandgrouse (family Pteroclidae) has alternatively been considered a member of Charadriiformes (shorebirds and allies [61, 64]) or Columbiformes [54, 57, 60]. Our reconstruction has the sandgrouse distantly related to both orders, and instead allied with Falconiformes. This relationship is unsupported elsewhere and we have little confidence in this placement. The novel placement of the tropicbird (family Phaethontidae) as sister to Sphenisciformes is similarly suspect.
Finally, we find no support in our mtDNA analyses for the neoavian clades 'Metaves' and 'Coronaves' identified from nuclear β-fibrinogen intron analyses , although our constraint tree allowed for this possibility (T Constraint; see additional file 1). A major difference between these trees involves the phylogenetic position of the perching birds (Passeriformes); while nuclear DNA analyses recover Passeriformes as a relatively derived clade within 'Coronaves' [57, 58], in T Optimal they instead comprise the basal-most lineage of Neoaves. This may be indicative of different phylogenetic signals in nuclear versus mtDNA sequences, as other mtDNA studies have also inferred a basal phylogenetic position of Passeriformes in Neoaves .
While T Optimal yields interesting hypotheses about avian relationships, the focus of this study is on estimating basal divergence times in Neornithes and we might regard topology as a nuisance parameter (and explicitly so in the BEAST T Flexible analyses). Topological error is usually not considered in divergence time estimation, but potentially could systematically bias age estimates through: (i) incorrect placement of fossil calibrations; and (ii) improper estimation of branch lengths. Through our joint consideration of T Consensus and T Optimal, we find that topology does have a systematic influence on inferred divergence times for nodes of interest (Table 3), but that for the present data set this influence differed in direction across methods and was generally insignificant when confidence/credibility intervals were considered. Removal of the restriction of a fixed topology in BEAST (T Flexible) yielded age estimates similar to those from Multidivtime analyses assuming T Optimal. Although yielding diffuse estimates, this 'relaxed topology' approach better reflects uncertainty in the underlying data.
An interesting result reported here is that rates of molecular evolution are found to be non-autocorrelated across the Neornithes tree (Tables 1 and 2), a result previously noted for virus and marsupial data sets . An autocorrelation of rates would imply an inheritance of the trait 'rate of evolution'. This makes intuitive sense when considering that ancestor and descendant lineages are likely similar in body size, generation time, DNA repair efficiency, population size and other traits influencing rates of sequence evolution, and the most popular molecular dating methods available indeed implicitly assume that rates are autocorrelated across a tree [66, 67]. However, even if 'rate of evolution' is heritable, nodes separated by long periods of time may accumulate sufficient rate variation that autocorrelation in this trait will break down. Alternatively, 'rate of evolution' may simply be more labile than we expect. Regardless, violation of an implicit autocorrelation assumption did not have significant effects on inferred dates for nodes of interest (Table 3). For example, r8s and Multidivtime, which each deal with rate variation in an ancestor-descendant fashion, deliver age estimates that match quite closely to those generated by BEAST, which does not make such an assumption.
All methods employed here agree that the basal divergences within Neornithes occurred in the Cretaceous (Table 3, nodes A-E), supporting the refutation of a Cenozoic origin of modern lineages [8, 9] mandated by the discovery of the 66 MY duck Vegavis iaai , which minimally forces five basal divergences into the Cretaceous. Moreover, our results are not dependent on this oldest fossil calibration, as analyses in r8s, PATHd8 and Multidivtime without using the Vegavis constraint returned nearly identical results to those reported here (data not shown); indeed, we must paradoxically conclude that this oldest undisputed neornithean fossil was essentially uninformative in our molecular dating analyses. Given the consensus across 'relaxed clock' methods employing very different assumptions about how molecular substitution rate evolves, we regard an Early Cretaceous origin of Neornithes as robustly supported. This inferred Cretaceous origin, and consequent survival of several avian lineages across the K-Pg boundary , is consistent with previous molecular studies [4, 16–26] and is supported by historical biogeography reconstructions .
An explanation occasionally offered for the antiquity of molecular dates is that rates of sequence change may speed up during radiations , consistent with a basic tenet of punctuated equilibrium theory where most character change is concomitant with speciation , possibly resulting from stochastic effects during founder effect speciation [52, 71]. If true, an elevated number of molecular substitutions recorded during diversification could erroneously be interpreted as a longer time span at a slower background rate, resulting in a systematic overestimation of divergence times for all nodes predating the radiation. Some evidence exists for a correlation between speciation and character evolution [52, 72–74], although a study of island radiations failed to detect such an effect . While punctuated molecular evolution may be less frequent in animals (18% of cases versus 44% and 60% in plants and fungi, respectively ), this effect is nevertheless a prime candidate to explain the strong and persistent discrepancy between molecular- and fossil-based divergence estimates. However, we find no signal for punctuated (speciational) molecular evolution  in the present data set. In addition, we fail to detect an accelerated rate associated with either the K-Pg boundary or during the initial diversification of Neornithes (Figure 5). If rates of sequence change were strongly influenced by diversification, we would expect clear departures from the inferred mean standardized substitution rate . Although Cenozoic rates tend to be more variable than Mesozoic (ancestral) rates, we find no evidence for an obvious acceleration in substitution rate associated with any of the major nodes for any genetic partition. Although these two approaches to identifying episodic evolution would ideally involve more comprehensive taxon sampling, if punctuated evolution were primarily responsible for inflating inferred molecular dates then we likely would have detected the effect with the present data matrix.
Rather than narrowing the formidable 'rock-clock' gap through application of methods designed to minimize biases and accommodate uncertainty, we have instead considerably reinforced it. In part, the discordant age estimates can be explained through reference to the genuine time-lag (δ True MA-FA ; see Figure 1) between the divergence of genetic lineages (predating speciation) and the evolution of diagnostic morphological characters (postdating speciation). However, given the estimated magnitude of δ Realized MA-FA (Table 3), it is unlikely that δ True MA-FA alone explains the dissonance. One the one hand, while the fossil record may be adequate for many evolutionary questions , there are clear instances where it is not. The 66 MY Vegavis iaai , for example, requires the coexistence of Paleognathae representatives; however, Cretaceous paleognaths are unknown. This may be a result of a geographical bias in fossil sampling favouring the northern hemisphere [2, 17, 69, 77–79], as many taxa are hypothesized as having southern hemisphere (Gondwana) origins . Although fossil gap analysis implies that a Cretaceous origin of several avian orders is unlikely , this method improperly assumes that fossils are randomly distributed (uniformly recovered through time) since the origin of a taxon; alternative fossil recovery curves can support very different scenarios, including scenarios that are more consistent with molecular genetic timelines , even when the fossil record is particularly sparse . Although rightly considered with caution, the increasing number of fragmentary remains of putative neornithean lineages from the Cretaceous  lends credence to the ancient origin and diversification of Neornithes. On the other hand, there may yet be some unrecognized biases in the analysis of molecular genetic sequences, and our results suggest new directions for future avian divergence time studies (described below).
Despite broad agreement on the antiquity of basal divergences in Neornithes, arbitration among macroevolutionary models [1–4] to best describe the history of more derived lineages is complicated by disparity amongst various 'relaxed clock' results. Two important points of distinction can be recognized (Figures 3 and 4). First, Dating5 (overdispersed clock) and PATHd8 (rate smoothing across sister lineages) both infer bursts of speciation across the avian tree, while r8s (smoothing in an ancestor-descendant fashion), Multidivtime (Bayesian modelling of ancestor-descendant autocorrelated rate evolution) and BEAST (Bayesian modelling of rate evolution without an autocorrelation assumption or fixed topology) instead infer a more gradual diversification of Neornithes. Second, PATHd8 alone supports a prominent radiation of avian families in the Cenozoic, a scenario statistically rejected in many instances by the remaining four analyses. Although published PATHd8 divergence time estimates for Neoaves using nuclear DNA produced similarly young estimates , a reanalysis of these same data using Multidivtime roundly refuted the findings , echoing the incongruence of PATHd8 reported here. While the better reconciliation between molecular and fossil age estimates that PATHd8 offers seems satisfying at first, the unique discrepancy of this method from the other much more rigorous and biologically realistic methods raises concern.
Given the apparent lack of autocorrelation of substitution rates, together with the intrinsic topological uncertainty in the Neornithes tree, the analyses in BEAST best reflect our current understanding of early avian evolution (Figure 4). Without the restrictive assumptions inherent in other 'relaxed clock' methods, these analyses unambiguously support a Cretaceous origin and diversification of basal avian lineages. Even when considering large inferred credibility intervals, 37 early avian divergences are restricted to the Cretaceous, similar to that found through the analysis of nuclear DNA sequences . It should be noted, however, that these nodes mostly represent order- and superfamily-level divergences; the majority of families sampled here (80 of 100 in BEAST) have their origins either overlapping or restricted to the Paleogene, consistent with interpretations from the fossil record . In this respect, our results are similar to those inferred from a comprehensive study of the tempo of early mammalian evolution . The results of both studies indicate that significant cladogenesis occurred in the Cretaceous, but that many crown groups diversified in the Cenozoic.
In regards to the 'rock-clock' debate [33, 49], we feel that much of the perceived dissonance between fossil- and molecular-based age estimates stems from a comparison of different phylogenetic entities: molecular phylogeneticists generally deal with stem-group origins, while paleontologists generally concentrate on crown-group diversification . Moreover, it is too rarely emphasized that when dating the same node a genuine discrepancy is expected, as coalescent times (T gene; see Figure 1) will predate cladogenesis while morphological diversification (T morphology; see Figure 1) will postdate cladogenesis. The reality is that in normal practice neither group directly addresses the main parameter of interest, the timing of speciation (T species; see Figure 1), the pattern of which is essential to the understanding of the origins of biodiversity. However, tools do exist to better estimate speciation times from both fossils (e.g. accounting for 'ghost' lineages ) and genetic data (e.g. explicitly modelling ancestral effective population sizes to account for coalescent times ). Further, molecular methods can be augmented with greater information from the fossil record, possibly by incorporating models of preservation bias into temporal constraints . Newly developed methods exist where probability distributions can be applied to calibrated nodes in a Bayesian framework [47, 84, 85]. Although we recognize that this approach is superior in that it can lend more credence to the fossil record than standard minimum-age constraints, we refrained from using such methods here as there is currently insufficient information with which to specify meaningful prior distributions for most avian diversification times. Realization of such distributions will require more communication between paleontologists and molecular phylogeneticists [86, 87].
One possible explanation for the discrepancy between molecular and fossil data in dating early divergences of avian lineages has been that the genetic data have been misinterpreted. In this vein, the ancient age estimates reported from previous molecular studies are seen as artefacts of oversimplified or improperly executed methods. Here we have examined this conjecture by accommodating uncertainty in genetic divergence time estimates. Through analyses of a large mtDNA matrix using multiple cladistically tested calibrations, alterative tree topologies and several sophisticated relaxed clock methods we have found that the molecular estimates are robust to varying assumptions about the evolution of evolutionary rates and consistent with those from previous studies. Our findings thus strongly support pre-K-Pg genetic origins for multiple modern bird lineages, including various extant orders and families, in contrast to the model of post-K-Pg diversification derived from a narrow interpretation of the fossil record.
Aligned fragment lengths of mtDNA sequences (total 4594 bp). Codon positions are combined across all protein-coding genes (ND1, ND2, cytochrome b), and RNA includes 12S and nine tRNA genes (L, I, Q, M, W, A, N, C, Y).
Aligned length (bp)
Variable sites (%)
Parsimony informative sites
First codon positions
Second codon positions
Third codon positions
Fossil calibrations. All internal calibrations for Neornithes are treated as minimum ages. External calibrations are treated as bounded (lower, upper) constraints. See  for fossil citations and details. See Figure 2 for placement of calibrations in the alternative trees.
Stem Upupidae + Phoeniculidae
Stem Coraciidae + Brachypteraciidae
Inferring dates on an incorrect tree topology will lead to estimates that are suspect, if not wholly incorrect. We seek to accommodate the existing uncertainty about avian phylogenetic relationships by dating nodes on two alternative trees. The first topology considered is taken from Figure 27.10 of Cracraft et al , which is a consensus tree based on previous molecular- and morphology-based phylogenetic reconstructions. This tree is conservative in that all represented branching events are well supported by independent lines of evidence; uncertain relationships among lineages are presented as hard polytomies. This topology is referred to as T Consensus. The second topology considered was derived from a partitioned-model maximum likelihood search using the program RAxML-VI-HPC 2.2.3 . The data were divided into four partitions: first, second and third codon positions of mitochondrial cytochrome-b, ND1 and ND2 genes, and concatenated mitochondrial 12S rDNA and tRNA genes (L, I, Q, M, W, A, N, C, Y). A collapsed consensus tree derived from Cracraft et al (thick branches only of Figure 27.10 in ) was used as a backbone constraint to make tree searches more efficient (T Constraint; see Additional file 1). Several hundred heuristic searches were performed under the GTRMIX model assuming a range of values for both the initial rearrangement setting (i; range 5–20) and number of rate categories (c; range 5–25). The topology of the maximum likelihood estimate (MLE) is referred to as T Optimal.
For the non-Bayesian dating methods, we accommodate uncertainty in branch length estimation through a non-parametric bootstrapping procedure. For each original data partition, 100 pseudoreplicate datasets were generated using the program SEQBOOT of the PHYLIP 3.6 package ; these bootstrapped matrices were concatenated to produce 100 pseudomatrices with the same size and number of partition-specific characters as the original matrix. For T Optimal, branch lengths and substitution model parameters were estimated using a partitioned GTR+G model in RAxML. For T Consensus, branch lengths and all substitution model parameters were estimated from each bootstrap matrix on the fixed topology using the GTR+I+G substitution model in PAUP*  because RAxML cannot evaluate a tree containing mutlifurcations. Using these procedures we generated 100 trees with branch lengths (phylograms) for each topology.
Owing to the concern that assumptions of particular analytical methods may systematically bias divergence time estimates, we compare several methods for accommodating among-lineage rate heterogeneity for the purpose of estimating the divergence times of modern avian lineages. To facilitate direct comparison, all methods utilize the same fossil calibrations and tree topologies outlined above.
First, the program r8s 1.7  was used to estimate rates and dates for internal nodes in the phylogeny via penalized likelihood (PL). This semiparametric procedure combines a parametric model that allows for different rates on each branch in the tree  with a non-parametric penalty function that penalizes rates that change too quickly along the tree from ancestor to descendent branches ; a smoothing parameter (λ) determines the relative contribution of the penalty function. The optimal smoothing value was determined individually for each topology through a sequence-based cross-validation procedure  using penalty = add and the normalized (χ2) errors. Confirmation of the optimal smoothing values was obtained via multiple optimizations with different initial conditions (by setting num_restarts = num_time_guesses = 3). Confidence intervals for node ages were determined using the distribution of estimated ages across bootstrapped trees assuming the optimal smoothing value from the original phylograms. Summary of the mean estimate and associated error for a given node was accomplished using the profile command in r8s.
Second, the program PATHd8 [98, 99] also smoothes rates across the tree, but does so by calculating an average path length from an internal node to all its descending terminals; deviations from a molecular clock are corrected through reference to fossil calibrations. Important distinctions from r8s above are that PATHd8 smoothes rates between sister groups (rather than from ancestor to descendent nodes) and it does this locally rather than over the entire tree. The same 100 phylograms as analyzed with r8s above were used to generate confidence intervals on divergence times, although summary statistics were calculated manually.
As a third approach, the Bayesian MULTIDISTRIBUTE package  represents an attempt to explicitly model rate change across a tree [67, 101, 102]. Here, the logarithm of the rate at the end of a branch is modelled with a normal distribution, the mean of which has an expected value equal to the rate at the beginning of the branch; thus, rates are implicitly assumed to be autocorrelated from ancestor to descendent nodes, although this autocorrelation may decay with increasing branch lengths. The posterior probability distribution of divergence times is approximated by samples from a Markov chain Monte Carlo (MCMC) chain. The data were partitioned as noted above. For each partition, estimates of the transition/transversion rate ratio and the gamma site class-specific rates under the F84+G model (the most complex model supported by the MULTIDISTRIBUTE package) were calculated in the baseml program of the PAML 3.15 package . The output from baseml was used as the input for the MULTIDISTRIBUTE program estbranches, which produced MLEs of branch lengths and their approximate variance-covariance matrix. Although differing in branching order, T Consensus and T Optimal had similar overall tree lengths; as a result, the same priors were applied to both topologies in the program Multidivtime: rtrate = rtratesd = 0.012, and brownmean = brownmeansd = 0.01. As one of our external calibration points bounds the root node (crocodile-bird split at 251-243 MY), date priors were less important and were defined narrowly (bigtime = 255 MY, rttm = 247 MY, rttmsd = 1.5 MY). The program was run without the assumption of correlated changes in rates across data partitions. Following a burnin of 106 samples, 104 samples were taken at a sampling interval of 102. All analyses were repeated with different inseed values to check for convergence of the MCMC chain.
Fourth, the overdispersed clock method of Cutler  represents a strong departure from other approaches to dating in the way it models branch length heterogeneity. Instead of treating a variable number of substitutions across lineages as indicative of changes in substitution rate across the tree, Cutler's method assumes that rate is stationary, but with high variance. Thus, rather than assuming that the number of substitutions along a lineage is Poisson distributed (where the variance is equal to the mean), the observed variation across branches is accommodated by the larger variance afforded through a Gaussian distribution. As a result, branches with either particularly high or low numbers of substitutions need not be clustered on the tree; in other words, heritability (phylogeny) plays no role in the observed numbers of substitutions. The program Dating5  calculates χ2 statistics for both the constant-rate Poisson and stationary models and compares the overall fit of the models using a likelihood ratio test. In addition, the program calculates R, the index of dispersion (the ratio of the variance to the mean number of substitutions) under the stationary model. Dating was restricted to T Optimal as the current version of Dating5 does not allow for polytomies. In addition, we could not achieve likelihood convergence with the partitioned data, and so reported results are from concatenated sequences. From asymptotic likelihood theory, 95% confidence intervals were calculated using a threshold of 2 log likelihood units around the MLE.
Finally, the program BEAST 1.4.6  differs in two important ways from the dating methods listed above. First, it does not require a fixed topology; rather, branch lengths, topology, substitution model parameters and dates can be estimated simultaneously. This incorporation of uncertainty in topology may be particularly important for the present data set, where relationships amongst many clades are uncertain. Second, BEAST does not assume that substitution rates are necessarily autocorrelated across the tree. Although intuitively satisfying, autocorrelation of rates has not been demonstrated in the literature; in fact, Drummond et al  report that many empirical data sets do not exhibit significant autocorrelation of rates.
BEAST 1.4.6 offers two statistical distributions for describing the change in rate across a branch : rates can be drawn independently from either an exponential distribution (UCED) or a lognormal distribution (UCLN). To find which distribution fit the present data best, we initially fixed the tree topology to T Consensus. The data were partitioned as above, with each partition assigned a GTR+I+G substitution model. BEAST MCMC runs of 25 × 106 generations following a burnin of 105 generations were performed for UCED, UCLN and CLOCK models. To arbitrate between models, we calculated Bayes factors by comparing harmonic mean model likelihoods . For both non-autocorrelated models, we also calculated the covariance among branch rates, which indicates the degree of ancestor-descendant autocorrelation of rates across the tree. Using the optimal model, we then accommodated topological uncertainty by removing the restriction of a fixed tree. However, we did constrain certain clades (the backbone constraints described above) to be monophyletic to facilitate the placement of calibration points and increase the efficiency of the MCMC search. Three replicate runs of 25 × 106 generations were performed to check for convergence of the MCMC chain. Mean parameter estimates and 95% highest posterior densities (HPDs) were determined through analyzing the combined BEAST tree files in TreeAnnotator 1.4.6 . We refer to these results with the topology T Flexible.
For each of the five dating methods above we calculated the average discrepancy between molecular (MA) and fossil (FA) estimated ages for those nodes that had fossil calibrations. The MA used in these calculations was the mean estimate from MCMC samples (Multidivtime, BEAST), or the optimal estimate from the empirical data matrix (r8s, PATHd8, Dating5). The value MA-FA is equivalent to δ Realized MA-FA described above, and gives an indication of the degree of agreement between the molecular data and the fossil record.
If present, episodic (or punctuated) molecular evolution could seriously bias molecular genetic estimates of divergence time by systematically overestimating the ages of all nodes that preceded the punctuation. We therefore tested for the presence of episodic evolution in two ways. First, we used the method of Pagel et al [52, 108] which quantifies the proportional contribution of punctuated (β) and gradual (g) evolution to path lengths in a phylogeny, based on extent of association between sequence change and cladogenesis events. This method requires a fully bifurcating tree, and so analyses were limited to our optimal reconstruction T Optimal. To test for this signature we analyzed maximum likelihood trees from RAxML bootstrap analyses (n = 100). Second, to test whether substitution rate accelerated coincident with or following the K-Pg boundary we plotted standardized inferred substitution rate (per data partition) versus inferred divergence time from the Multidivtime analyses above. If the K-Pg boundary extinctions facilitated diversification of avian higher-level taxa, it could produce a spike in this plot  detected as a departure from the mean standardized rate. These two tests are complementary in that the first focuses specifically on the effects of speciation, whereas the second focuses on absolute time effects, which may or may not be related to speciation.
As in the case of episodic evolution, a simple lack of historical signal present in molecular sequences could generate erroneous divergence time estimates. We therefore assessed the 'information content' present in our mtDNA matrix by regressing posterior 95% credibility width against posterior mean divergence time. When the when the amount of data is saturated then this regression will produce a linear relationship (i.e. R 2 → 1), the slope of which indicates the amount of uncertainty that cannot be reduced through adding longer sequences [7, 85], but can be reduced through adding independent loci.
corrected Akaike Information Criterion
strict molecular clock
highest posterior density
maximum clade credibility
Markov chain Monte Carlo
maximum likelihood estimate
time to the most recent common ancestor
uncorrelated exponential distribution
uncorrelated lognormal distribution
We thank A Stamatakis (RAxML), J Thorne (MULTIDISTRIBUTE) and D Cutler (Dating5) for assistance with their respective software packages, and S Ho for assistance with the new methods available in BEAST. M van Tuinen offered indispensable advice regarding an initial set of fossil constraints and RB Payne provided essential comments on an earlier (encyclopaedic) draft of this manuscript. Three anonymous reviewers offered critical suggestions to improve the manuscript. JWB thanks R Pollard for sustained encouragement throughout. We thank A Baker and A Lindsay for assistance with laboratory work. Funding was provided by the National Science Foundation and the University of Michigan.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.