|Home | About | Journals | Submit | Contact Us | Français|
The Early Cretaceous is a critical interval in the early history of birds. Exceptional fossils indicate that important evolutionary novelties such as a pygostyle and a keeled sternum had already arisen in Early Cretaceous taxa, bridging much of the morphological gap between Archaeopteryx and crown birds. However, detailed features of basal bird evolution remain obscure because of both the small sample of fossil taxa previously considered and a lack of quantitative studies assessing rates of morphological evolution. Here we apply a recently available phylogenetic method and associated sensitivity tests to a large data matrix of morphological characters to quantify rates of morphological evolution in Early Cretaceous birds. Our results reveal that although rates were highly heterogeneous between different Early Cretaceous avian lineages, consistent patterns of significantly high or low rates were harder to pinpoint. Nevertheless, evidence for accelerated evolutionary rates is strongest at the point when Ornithuromorpha (the clade comprises all extant birds and descendants from their most recent common ancestors) split from Enantiornithes (a diverse clade that went extinct at the end-Cretaceous), consistent with the hypothesis that this key split opened up new niches and ultimately led to greater diversity for these two dominant clades of Mesozoic birds.
The recent availability of larger phylogenetic datasets and numerical approaches that account for temporal variation in sampling have led to an explosion in quantitative approaches to exploring macroevolution in palaeontology [1–5]. Among these the dinosaur–bird transition has seen considerable attention. For example, recent comparative studies indicate that miniaturization along the lineage leading to birds was vital to their Mesozoic radiation and that birds evolved faster than other non-avian dinosaurs [4–7]. However, these studies focused on single continuous measures (limb proportions and body size), or contained limited numbers of Mesozoic avian taxa, insufficient for exploring detailed features of early bird evolution. Over the past three decades, our understanding of early avian evolution has been greatly advanced by discoveries from Mesozoic localities, particularly the Lower Cretaceous Jehol Biota [8–10]. These exquisitely preserved Early Cretaceous bird fossils have narrowed the large morphological gap between the iconic ‘Urvogel’ Archaeopteryx and living birds, and importantly indicate that many morphological features once considered unique to crown birds had already evolved in stem-group taxa only slightly younger than Archaeopteryx [8,11–13]. For instance, the long-bony tail characteristic of non-avian theropods (plesiomorphically retained in Archaeopteryx and Jeholornis) has been abbreviated and replaced by a pygostyle in Confuciusornithiformes and Sapeornithiformes [14–16]; derived features present in crown-group birds, e.g. the keeled sternum, fused compound bones, alula, and fan-shaped rectrices, are well established in Enantiornithes and Ornithuromorpha, which are the most successful clades of Mesozoic birds [12,17–20]. Enantiornithes and Ornithuromorpha can be easily distinguished from earlier diverging (basal) birds, but morphological disparities among their own Early Cretaceous members are relatively small [21,22]; therefore, it is likely that a burst of phenotypic innovation took place close to the origin of Enantiornithes and Ornithuromorpha, and facilitated their adaptive radiation. Such rapid morphological evolution is apparently followed by a deceleration in their subsequent Early Cretaceous history. However, the quantitative analysis of morphological tempo between primitive avian lineages has not previously been performed. Here we use phylogenetic methods and sensitivity tests to discern patterns of morphological evolutionary tempo within Early Cretaceous avian phylogeny for the first time.
Phylogenetic analysis was carried out using the recently published comprehensive morphological characters of Mesozoic birds in Wang et al. . Tree searches were performed in PAUP (v. 4.0b10, ), using the parsimony method and the settings of Wang et al. . Four most parsimonious trees (MPTs) were recovered (electronic supplementary material, figure S1), a small enough number that all trees were subjected to subsequent rate tests to gauge phylogenetic uncertainty. The complete character description, matrix, and MPTs are provided in the electronic supplementary material.
Before rate analyses could proceed it was necessary to time-scale our trees as branch durations (along with completeness, see below) form the denominator in our rate equation. This first required establishing ages for terminal tips (taxa). However, the absolute duration of a Mesozoic avian species is rarely known with certainty, and in most cases only the geological age of the sediments from which the taxon was collected is known [24,25]. Therefore, the duration of the fossil species used here represents stratigraphic uncertainty rather than a true range. Ages for our taxa (given in the electronic supplementary material) were collected from the primary literature and are expressed as lower (first appearance datum, FAD) and upper (last appearance datum, LAD) bounds of the geologic stage(s) or epoch(s) to which the fossil-containing sediments can be referred. However, direct dates such as radioisotopic dating, were used where available, e.g. for the Huajiying, Yixian, and Jiufotang formations [26,27].
Traditionally, trees of fossil-only tips have been scaled using a ‘conservative’ approach where each internal node shares the same date as its oldest descendant, producing multiple zero-length branches, ZLBs [28,29]. This approach prevents rate calculation as it forces many denominators to be zero, effectively implying infinite rates. ZLBs may be removed by adding some constant, for example, applying the ‘minimum branch length’ (mbl) method in Laurin , which forces all branches to have a fixed minimum duration (here we used 1 Ma). Alternatively, the ‘equal’ method may be used, which iteratively removes ZLBs by allowing them an equal share of the duration of the first preceding (ancestral) positive length branch . Here we applied both methods as a sensitivity test. The ‘mbl’ method was applied to the strict consensus tree (for visualization; figure 1) and each MPT (for rate analysis; figure 2) using the timePaleoPhy function in the R package paleotree . Similarly, the ‘equal’ method was applied to each MPT using the DatePhylo function in the R package strap , using a ‘root-length’ (required for the base of the tree where there is no preceding positive length branch) of 2 Ma. In addition, a date randomization approach , comparable with that in , was also applied to incorporate the uncertainty involved in dating terminal taxa (i.e. the gap between an FAD and an LAD were correctly treated as uncertainties rather than ranges). Here we performed 100 replicates—a value that served as a compromise between computation time and a full delineation of dating uncertainty—of dating randomization, where taxa were assigned a numerical date drawn randomly from a uniform distribution bounded by its FAD and LAD.
In addition to the above sensitivity tests (four different MPTs, two different time-scaling methods, 100 different dating randomizations), we produced trees with Late Cretaceous taxa both included and excluded. In total, this amounts to 1 600 separate time-scaled phylogenetic hypotheses for Mesozoic birds, all of which were subjected to the same rate analysis protocol described below.
At its simplest our evolutionary rate calculation is the number of changes (along a branch or within a clade) over the duration of that branch (or the sum of all branches within a clade). Thus, high rates can be found for very low numbers of character changes (due to a short duration branch or branches), and conversely, low rates may be found where there are many character changes due to a long duration branch (or branches). The importance of branch duration is why we subjected our data to the multiple time-scaling approaches described above.
However, our rate equation really has two terms in its denominator, time and completeness. Because our morphological data are morphological characters, and character states frequently cannot be coded for all taxa (either due to incompleteness or inapplicability) accounting for missing data is of key importance. In practice, a branch's completeness corresponds to the number of characters for which states can be scored at both ends, and hence the number of characters for which a change could be observed . A highly complete branch offers multiple opportunities to observe a change, and a highly incomplete branch few opportunities. Thus, without accounting for this difference we may, incorrectly, assign a higher rate to a more complete branch and a lower rate to a less complete branch even if the true pattern is homogeneous rates. Here, missing data are dealt with in two ways. Firstly, completeness is added to the denominator of our rate equation to account for this issue. Thus, per-branch and per-clade rates explicitly include completeness and are measured as changes per character per lineage million years. However, such raw rates should not (and here, are not) read as direct evidence for high or low rates. This is because when completeness is particularly low extreme values can become common, but we should not have confidence in their reality. As an example, if only one binary character can be observed along a branch then we either observe one change (the maximum possible number) or no changes (the minimum possible number). This can lead to an apparently very high rate in the first case, or an absolute minimum (zero) rate in the second case. However, in reality, a single character should give us no confidence in determining the significance of the high or low rate observed. Secondly, then, completeness is used in our statistical tests effectively as our sample size and helps avoid misinterpreting absolute rates.
In practice, our rate calculations require the input of three separate trees , with branch lengths reflecting: (i) number of character changes, (ii) completeness (number of opportunities to observe changes), and (iii) time. At this stage, we have only generated the trees with branch lengths relating to time. In an earlier protocol , the number of character changes was determined by outputting data directly from PAUP , using a parsimony algorithm to generate ancestral states, then simple counts of differences between the states at either end of a branch were performed. However, this was problematic as the algorithms used (accelerated transformation (ACCTRAN) and delayed transformation DELTRAN)) did not adequately deal with missing data, filling in values for every character at every internal node even for clades where all tips exhibited the missing state (?). This biased incomplete branches to appear to confidently exhibit low rates, when in reality confidence should have been lower and rates higher. It did, however, simplify the generation of the completeness tree, with all internal branches being considered 100% complete and all terminal branches as complete as the tip to which they lead. Another problem with filling missing states is that this is illogical for inapplicable characters, for example, estimating feather colour in a featherless taxon.
Ancestral states can better be captured by likelihood methods that employ branch lengths in their calculation [33,34], and offer more realistic models of evolution [35–37]. Such maximum-likelihood methods consider branch lengths (reflecting time) in ancestral state reconstruction in such a way as to make the character states observed in the terminal taxa most probable [33,38]. Here, we specifically applied the method of Yang et al. , as implemented in the rerootingMethod function in the R package phytools . Once the ancestral states are reconstructed, using the above function where data are available (and collapsing ambiguities into polytomies) or simply placing question marks (missing values) where it is not, both the number of character changes along each branch and the completeness of each branch can be calculated.
Because of the issue regarding completeness, we followed previous studies [5,32] in only asking whether rates are significantly high or low, rather than what their absolute values are (although these are shown in the electronic supplementary material). More specifically we ask, for every branch or clade in the tree, whether that branch or clade has a higher or lower rate of morphological evolution than the rest of the tree. We do this using the likelihood tests outlined in Lloyd et al.  and refer the interested reader there for a more complete description. In short, the algorithm models character change using a Poisson process (all changes are treated as independent) and then applies a likelihood ratio test (LRT) between two models, a single rate model (one rate occurs across the entire tree) versus a two rate model (the rate on the target branch, or within the target clade, is different to that for the rest of the tree). For each individual test (i.e. a specific branch or clade), the LRT returns a p-value that can then be compared with a user-defined α (here we use 0.01) to decide whether the result is significant (i.e. if the p-value falls below our α then the two rate model is preferred). Owing to the large number of tests being performed (every branch or every clade), a multiple comparison correction is required to avoid type I errors (false positives) and here the algorithm uses the Benjamini and Hochberg false discovery rate test .
An important additional bias to consider here is the difference between terminal and internal branches. Specifically, using current time-scaling methods (outlined above) internal branches are often found to be shorter than terminal branches, and consequently, tend to exhibit higher rates . It is not yet clear why this is the case, or whether it might be a more general problem. However, in order to avoid results driven purely by this difference all significance calculations were repeated for either internal- or terminal-only branches. In other words, for an internal branch its rate was compared with only other internal branches, or for a clade either only the terminal or only the internal branches were included in its rate calculation and subsequent comparison. Note that in the latter case ‘cherries' (clades containing only two terminal taxa) could not be tested as by definition they contain no internal branches.
As a final sensitivity test, we considered the inclusion or exclusion of Late Cretaceous fossil birds. These are known from fragmentary materials and their affinities have never been comprehensively explored in the context of a phylogenetic analysis. This can be demonstrated quantitatively by comparing the mean character completeness of the better known Late Cretaceous birds included here (43%) with that of our broader Early Cretaceous sample (50%; electronic supplementary material figure S2). We therefore restrict our main analyses to the Early Cretaceous by using the timeSliceTree function in the R package paleotree . By applying a cut-off of 113 Ma (the top of the Albian), we effectively focus on the first 39 Ma of avian evolution (152 Ma–113 Ma). However, the ‘removed’ younger taxa were still considered in both the tree searches and time-scaling. They were also retained for a separate set of evolutionary rate tests to gauge whether their exclusion really affects our results.
The time-scaled phylogenetic tree suggests that the divergence between basal avian clades more derived than Archaeopteryx occurred at the Jurassic–Cretaceous boundary (figure 1), using the ‘mbl’ scaling method . Applying the ‘equal’ method  instead pushes the same nodes back into the Late Jurassic, e.g. Pygostylia (150.4 Ma) and Ornithothoraces (145.5 Ma; electronic supplementary material, figure S3), much earlier than the oldest known record of either clade—the late Early Cretaceous (130.7 Ma) [17,19,41]. Thus, in this case, the ‘mbl’ method is more conservative, minimizing the amount of inferred missing history.
For all 1 600 trees examined, a null test for homogeneous rates across the whole phylogeny was confidently rejected (p < 10−43 for the ‘mbl’ trees and p < 10−54 for the ‘equal’ trees). Thus, individual branch and clade tests were warranted and results for the first MPT are shown in figure 2 (‘mbl’; first MPT) and figure 3 (‘equal’; first MPT), and full results on Dryad. Including Late Cretaceous birds necessarily lowers absolute rates, because it requires the inclusion of longer branches (figure 1). However, significance tests for branches and clades in the Early Cretaceous are not substantially altered by the inclusion of Late Cretaceous taxa (cf. figure 2, electronic supplementary material, figure S4, and figure 3, electronic supplementary material, figure S5; see Dryad for complete results), showing that our results are not biased by excluding the Late Cretaceous lineages.
Using both time-scaling methods, branches subtending primitive taxa, including the long-bony tailed Archaeopteryx, Jeholornis and the basalmost pygostylians Confuciusornithidae, consistently have non-significant rates (neither significantly higher nor lower than the average rate across the whole phylogeny) whether branches are treated together or separately as internals and terminals (figures 2 and and3;3; Dryad data). The branch subtending Sapeornis, the direct outgroup to Ornithothoraces, has significantly low rates where branches are treated together, but is non-significant when treated separately (figures 2 and and3;3; Dryad data).
When branches are compared using the ‘mbl’ method, significantly high rates dominate the branches subtending Ornithothoraces, Enantiornithes, Ornithuromorpha, and certain branches within these latter two groups (figure 2a). However, these high rates are severely diminished when internal and terminal branches are treated separately, with the majority of the branches having non-significant rates, in particular the branches subtending Enantiornithes and Ornithuromorpha (figure 2b). However, certain internal and terminal branches still retain significantly high rates over half of the dating replicates: the internal branch subtending Pengornithidae, and the terminal branches subtending Protopteryx, Piscivoravis, and Yanornis (figure 2b). Over half of the branches within the Hongshanornithidae have non-significant rates and the terminal branches subtending Hongshanornis, Parahongshanornis, and Tianyuornis tend to be dominated by significantly low rates (figure 2).
Branch likelihood tests using the ‘equal’ method produce essentially the same results, except for minor differences regarding the branch subtending Ornithothoraces (figure 3; Dryad data). Significantly high rates are found along the branch subtending Ornithothoraces using both the ‘mbl’ and ‘equal’ methods when branches are treated together (figure 3). Overall, this suggests a high rate of morphological change when Ornithothoraces diverged from Sapeornithiformes.
The clade likelihood tests indicate that high rates dominated early diverging clades, with later diverging clades having non-significant rates when using the ‘mbl’ time-scaling method (electronic supplementary material, figures S6 and S7; Dryad data); significantly high rates are more frequently observed using the ‘equal’ method (electronic supplementary material, figure S7; Dryad data). Within Enantiornithes, non-significant rates are commonly observed, and only certain derived clades have significantly high rates when the ‘mbl’ method is applied; whereas, approximately half of enantiornithine clades have significantly high rates using the ‘equal’ method. Pronounced rate variability is present within Ornithuromorpha using both the ‘equal’ and ‘mbl’ time-scaling methods: significantly low rates are discovered in all dating replicates within the Hongshanornithidae; by contrast, the sister clade to the Hongshanornithidae has significantly high rates over half of the dating replicates (electronic supplementary material, figures S6 and S7); overall, the Hongshanornithidae represent the earliest avian clade to consistently exhibit significantly slow rates. Enantiornithes are the most speciose group of Early Cretaceous birds , but significantly high rates are not widely distributed across this group: about half of enantiornithine clades have significantly high rates using the ‘equal’ method, and less than half using the ‘mbl’ method, which may explain why Early Cretaceous enantiornithines are generally uniform in morphology.
Ornithuromorphs and enantiornithines were more taxonomically diverse than their Early Cretaceous contemporaries, e.g. Sapeornithiformes and Confuciusornithiformes [22,43], but how this difference arose is unknown. Several newly reported ornithuromorphs from the Jiufotang Formation, including Archaeorhynchus, Schizooura, and Bellulornis, have begun to close the morphological gap between Enantiornithes and Ornithuromorpha [44–46], suggesting rapid rates of speciation and morphological diversification during the course of their early history. This is supported by the presence of significantly high rates found at the base of Ornithothoraces (figures 2 and and3).3). In particular, significantly high rates along the branch subtending Ornithothoraces and the non-significant rates along the early diverging branches within Ornithuromorpha and Enantiornithes warrant further discussion. Such differences can be explained by (i) an early burst of morphological evolution consistent with an adaptive radiation among early Ornithothoraces; or alternatively (ii) is an artefact of the time-scaling method used (the duration of the branches subtending Ornithuromorpha and Enantiornithes were underestimated and caused inflated rates). We compared the results of our branch likelihood tests using both the ‘mbl’ and ‘equal’ methods to test these two explanations. The ‘mbl’ time-scaling method sets the ZLBs to a fixed minimum duration of (here) 1 Ma, and the ‘equal’ method iteratively removes ZLBs by allowing them to share duration equally with preceding non-zero length branches, making its effects less predictable. In this case ZLBs ‘smoothed’ by the ‘mbl’ method were shorter than those ‘smoothed’ by the ‘equal’ method (at least 50% shorter in our data). This resulted in higher absolute rates and a greater number of significantly high rates for the ‘mbl’ method. However, a simple prediction of lower rates across the tree for the ‘equal’ method is not borne out. Our branch likelihood test results show that significantly high rates (about half of the dating replicates) occur along early diverging branches within Ornithuromorpha when using the ‘equal’ method (figure 3 and Dryad data). However, instead of finding higher absolute rates, more significantly high rates, or a greater frequency (across dating replicates) of high rates using the ‘mbl’ method, non-significant rates were more common along the same branches (figure 2). Branch likelihood tests using both the ‘mbl’ and ‘equal’ methods produce essentially the same result for Enantiornithes. Therefore, the heterogeneous evolutionary rates between the branches subtending Ornithothoraces and early diverging branches within Ornithuromorpha and Enantiornithes are not likely to be biased by underestimation of branch length, and here the first explanation is considered more plausible. The transition from early rapid change to subsequent slower rates in the early history of Enantiornithes and Ornithuromorpha suggests a potential adaptive radiation (early burst) [47,48]. Morphological characters adaptive to arboreal (enantiornithines) and lake-shore habits (ornithuromorphs) facilitated early diverging members of both clades to invade new niches and thus further drive their morphological diversification [17,49]; however, as ecological niches became increasingly saturated, the rate of morphological change declined [50,51].
Recent comparative phylogenetic studies on body size and limb proportion changes across the dinosaur–bird transition reveal that reduced body size and the acquisition of novel features pertaining to flight have contributed to the sustained morphological and ecological radiation of birds [4,6,7,52]. Although these results came from large phylogenies containing many non-avian dinosaurs, significant shifts in limb proportions and body size in both Enantiornithes and Ornithuromorpha are suggested. In addition, these changes facilitated the occurrence of novel morphological features [4,7], which partially supports our result that high evolutionary rates occur at the base of Ornithothoraces. We divided the morphological characters into six anatomical subregions—skull, axial skeleton, pectoral girdle, forelimb, pelvis, and hindlimb, and then compared character changes during several major transitions, i.e. branches subtending Aves, Pygostylia, Ornithothoraces, Enantiornithes, and Ornithuromorpha. The result shows different patterns of morphological change during these transitions (figure 4). For instance, character changes relating to axial skeleton and forelimb peak along the branches subtending Aves and Pygostylia, whereas the pectoral girdle is more conservative in morphology during these transitions. When Ornithothoraces diverged from Sapeornithiformes, more than half of the morphological changes inferred occur in the pectoral girdle and forelimb, contributing to the refinement of their flight capability. Substantially different patterns of morphological modifications are observed between Enantiornithes and Ornithuromorpha: enantiornithines show numerous changes in their pectoral girdle and forelimb morphology, but their axial skeleton and pelvis are relatively conservative. By contrast, only a few character changes took place along the branch subtending Ornithuromorpha, and those modifications were more evenly distributed in those anatomical subregions. Collectively, our results shows that the pectoral girdle and forelimb morphology changes more than any other body parts, suggesting that aerodynamics played an important role in driving early avian evolution.
Overall, accelerated evolutionary rates occurred at the point when Ornithuromorpha split from Enantiornithes, in company with many synapomorphies evolved in each clade, e.g. elongated minor metacarpal, Y-shaped furcula, and straight scapula in Enantiornithes; keeled sternum, expanded proximal phalanx of major digit, modern-shaped pygostyle, and triosseal canal in Ornithuromorpha [8,12,17,18,53]. Such rapid evolutionary rates may have built an adaptive advantage that opened up new niches and ultimately led to greater diversity. This conclusion is supported both by our results and other studies, which indicate that Enantiornithes and Ornithuromorpha evolved rapidly to fill completely different ecological niches: enantiornithines are arboreal and ornithuromorphs are terrestrial [49,54–57]. In addition, sexual selection probably contributed considerably to the greater diversity of Ornithothoraces, given the complex feather morphologies present in this clade [58,59]. However, more data are needed, in particular for species younger than the Albian, to enrich our understanding of evolutionary patterns in Mesozoic birds.
Electronic supplementary material is available online, and additional data are available in the Dryad repository: http://dx.doi.org/10.5061/dryad.c128h.
M.W. designed the study; G.T.L. wrote scripts in R; M.W and G.T.L. ran phylogenetic and rates analyses; M.W. and G.T.L. wrote the manuscript.
The authors declare no competing interesting.
M.W. was funded by the National Natural Science Foundation of China (41502002), the Youth Innovation Promotion Association (CAS, 2016073). G.T.L. was funded by Australian Research Council grant no. (DE140101879).