Whether dinosaurs were in decline before their final extinction 66 Mya has been debated for decades with no clear resolution. This dispute has not been resolved because of inappropriate data and methods. Here, for the first time to our knowledge, we apply a statistical approach that models changes in speciation and extinction through time. We find overwhelming support for a long-term decline across all dinosaurs and within all three major dinosaur groups. Our results highlight that dinosaurs showed a marked reduction in their ability to replace extinct species with new ones, making them vulnerable to extinction and unable to respond quickly to and recover from the final catastrophic event 66 Mya.

Whether dinosaurs were in a long-term decline or whether they were reigning strong right up to their final disappearance at the Cretaceous–Paleogene (K-Pg) mass extinction event 66 Mya has been debated for decades with no clear resolution. The dispute has continued unresolved because of a lack of statistical rigor and appropriate evolutionary framework. Here, for the first time to our knowledge, we apply a Bayesian phylogenetic approach to model the evolutionary dynamics of speciation and extinction through time in Mesozoic dinosaurs, properly taking account of previously ignored statistical violations. We find overwhelming support for a long-term decline across all dinosaurs and within all three dinosaurian subclades (Ornithischia, Sauropodomorpha, and Theropoda), where speciation rate slowed down through time and was ultimately exceeded by extinction rate tens of millions of years before the K-Pg boundary. The only exceptions to this general pattern are the morphologically specialized herbivores, the Hadrosauriformes and Ceratopsidae, which show rapid species proliferations throughout the Late Cretaceous instead. Our results highlight that, despite some heterogeneity in speciation dynamics, dinosaurs showed a marked reduction in their ability to replace extinct species with new ones, making them vulnerable to extinction and unable to respond quickly to and recover from the final catastrophic event.

Nonavian dinosaurs met their demise suddenly, coincident with the Chicxulub impact in Mexico around 66 Mya; however, whether there was any long-term trend toward declining diversity leading to the Cretaceous–Paleogene (K-Pg) boundary has been controversial and debated for decades (1⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–14). This longstanding dispute has been prolonged partly because of differences in fossil datasets from different parts of the world and difficulties in rock dating but most importantly, methodological weaknesses—previous attempts have been nonphylogenetic, and analyses were conducted on simple time-binned tabulated data, resulting in a lack of statistical rigor (phylogenetic and temporal nonindependence have not been considered), and did not truly investigate evolutionary dynamics, such as speciation and extinction rates. In fact, patterns of speciation and extinction in dinosaurs have gone largely unstudied (8). Here, we study speciation dynamics (relationship between speciation and extinction rates) using an exclusively phylogenetic approach in a Bayesian framework.

If speciation and extinction rate were constant (but speciation was higher), we would expect to see a linear increase through time in the logarithm of the number of speciation events along each path of a phylogenetic tree (linear) (Materials and Methods and Fig. 1, model A). If speciation rate decreased through time but remained above extinction rate, then we would expect a curvilinear relationship (Fig. 1, models B and C). Such a relationship would reach an asymptote (speciation equals extinction) (Fig. 1, model B) and eventually, turn down as extinction rate surpasses speciation during the evolutionary history of the clade (Fig. 1, model C). The latter would correspond to a long-term pre–K-Pg demise in the case of dinosaurs. The distinction between such evolutionary dynamics can only be made using phylogenies with taxa sampled through time.

Theoretical models of speciation through time. If speciation and extinction rate were constant through time (but speciation was higher) in dinosaurian history, we would expect to see a linear increase through time in the logarithm of the number of speciation events along each path of a phylogenetic tree (model A). If speciation rate decreased through time but remained above extinction rate, then we would expect a curvilinear relationship (models B and C). Such a relationship would reach an asymptote (speciation equals extinction; model B) and eventually, turn down as extinction rate surpassed speciation during the evolutionary history of the clade (model C). The latter would correspond to a long-term pre–K-Pg demise in the case of dinosaurs.

Our results show that dinosaurs were in decline for a much longer period than previously thought—extinction rate surpassed speciation rate at least 40 My before their final extinction. This prolonged demise leaves plenty of time for other animal groups to radiate and flourish as more and more ecological niches open up, most prominently the pre–K-Pg expansion of crown mammals ( 46 ). Although Mesozoic dinosaurs undoubtedly dominated the terrestrial megafauna until the end of the Cretaceous, they did see a reduction in their capacity to replace extinct species with new ones, making them more susceptible to sudden and catastrophic environmental changes, like those associated with the asteroid impact.

Our study represents the first, to our knowledge, explicitly phylogenetic statistical treatment of speciation dynamics in dinosaurs. Unlike previous nonphylogenetic attempts to study changes in dinosaur taxic diversity across geological time bins ( 8 ⇓ – 10 , 13 , 35 , 44 , 45 ), our method is robust to sampling and other potentially confounding factors ( SI Text and Tables S1–S10 ) and can statistically detect decreases in net speciation, which is difficult if not impossible to establish using conventional methods. Furthermore, by accounting for the effects of shared ancestry, we provide a more accurate picture of dinosaurian speciation dynamics than the simple summing of species records through time.

Although we cannot positively identify a causal mechanism for the speciation downturn in dinosaurs, there are many possible global phenomena that occurred during the Cretaceous Period [e.g., the continued breakup of the supercontinents Laurasia and Gondwana (limiting free movement and eventual para- or peripatric speciation), intense prolonged volcanism ( 36 ), climate change ( 37 ⇓ – 39 ), fluctuations in sea levels ( 34 , 40 ), and ecological interaction with rapidly expanding clades ( 41 )]. To accurately identify causal mechanisms of Mesozoic dinosaurian demise, we recommend that future studies focus on a longer time period than just the last 10–20 My of the Cretaceous ( 4 , 13 , 42 , 43 ). In addition, our results highlight the importance of considering the expected increase in species number as clades expand and accounting for shared ancestry using phylogenetic approaches.

We can indirectly assess the influence of geography, such as segregation by geographic barriers ( 30 ), using Mesozoic eustatic sea-level reconstructions ( 34 ) as an additional covariate in our models. Although various hypotheses have been proposed regarding the influence of sea level on biodiversity in dinosaurs ( 35 ), the most compelling suggests that increasing sea level results in fragmentation of large landmasses and can alter geographical distributions of habitats. In turn, continental fragmentation can lead to geographical segregation, reproductive isolation, and ultimately, speciation ( 30 ). Our results, for the first time to our knowledge, support this hypothesis—we find a significant positive effect of sea level on speciation [ΔDIC (five-group quadratic − five-group + sea-level models) >16; the proportion of the posterior distribution of the Markov-chain Monte Carlo (MCMC) estimates that crosses over zero multiplied by two (p-MCMC) <0.0010] ( Tables S1–S10 )—although the effect is small; for every 1-meter increase in sea level, speciation events increased by 0.2–0.25%. Horner et al. ( 29 ) observed that the emergence of transitional morphotypes coincides with marine transgressions in Late Cretaceous rocks of western North America, consistent with our finding that rising sea levels induce speciation. Importantly, the inclusion of sea level in any of our models does not diminish the temporal decline in species proliferation, despite the substantial rise of sea levels worldwide by some 150–200 m throughout the Cretaceous ( Tables S2–S10 ).

An ecological limit on speciation or the filling of available niches ( 30 , 31 ) is commonly invoked to explain speciation slowdowns. Members of the same clade are more likely to compete for similar if not the same ecological niche or portions of ecospace ( 32 , 33 ), and the more numerous the number of contemporary lineages, the fewer the number of available niches. We tested such an effect by including a measure of intraclade niche competition—cladewise lineage diversity or the number of contemporary branches (including internal branches) for each taxon—in the model ( SI Text ). However, we find that cladewise lineage diversity is not significantly associated with speciation, and it does not explain the observed downturn; physical restrictions, such as geography or range sizes, could be more important.

Our results showing high levels of speciation in hadrosauriforms and ceratopsids, although consistent with previous findings ( 8 ), seem to contradict more recent work that suggests that these groups underwent a decline in morphological diversity during the last two stages of the Cretaceous of North America ( 13 , 26 ). These dinosaur species are morphologically and ecologically (at least at the family level) conserved ( 27 ), with most of the derived characteristics concentrated in their crania ( 24 ). Speciation can be high in these groups, despite the potentially low morphological diversity, because Cretaceous dinosaurs exhibited increased provincialism ( 28 ) (speciation arising from geographic isolation rather than sympatric niche partitioning), increased α-diversity (many more species with subtly varying skulls but identical postcrania sharing the herbivorous ecospace in single localities), and changing taxonomic composition of stable ecological community structures [ecological niches remain constant, but taxa filling those niches changed through time ( 27 , 29 )].

Ornithischians show a similar increase to theropods (∼0.06 speciation events for every 1 My) to ∼192 Mya followed by a slowdown to ∼114 Mya, at which point extinction rate exceeds speciation rate ( Figs. 2B and 3 ). Key morphofunctional features in oral food processing distinguish hadrosauriforms and ceratopsids from other ornithischians, permitting them to exploit major new food sources ( 23 , 24 ). Whether these herbivores were exploiting the new, small, fast-growing herbaceous angiosperms that became common and widespread as early as the Aptian–Albian (125–109 Mya) of the Early Cretaceous ( 25 ) is much debated. The powerful jaws and massive dental batteries of these herbivores might have been adapted to other tougher nonangiosperm plant food, and they benefited from a new adaptive complex in food processing.

Because birds underwent a radiation in the Early Cretaceous after their appearance in the Middle to Late Jurassic, one might expect that their pattern of speciation would be distinct from that of nonavian theropods. However, when we allow separate coefficients (intercept, slope, and quadratic terms) in our model to be estimated for birds and nonavian theropods, the resulting regression parameters were not significant: in other words, the speciation dynamics in Mesozoic birds are not distinct from those of nonavian theropods ( Tables S1–S10 ). This result is in line with recent findings of a high, sustained rate of change from the Late Triassic to the Early Cretaceous in the entire theropod lineage leading to Archaeopteryx and among the earliest birds ( 21 , 22 ).

Speciation in theropods follows a slower increase (∼0.07 speciation events for every 1 My) with an early onset of speciation slowdown from the Late Triassic ∼215 Mya to the Early Cretaceous ∼120 Mya, when extinction rate exceeds speciation rate ( Figs. 2B and 3 ). Although Theropoda contains one of the most morphologically diverse dinosaurian clades, the coelurosaurs, which include the giant carnivorous tyrannosaurs, parrot-like oviraptorosaurs, large pot-bellied therizinosaurs, ostrich-like ornithomimosaurs, small sickle-clawed dromaeosaurs, and birds (most of which are Cretaceous in age), they originated in the Early to Middle Jurassic ( Fig. 3 ), much earlier than expected from apparent fossil occurrences ( 8 ). Clades appearing even earlier (e.g., ceratosaurs, megalosauroids, and allosauroids) also persist into the Late Cretaceous, all of which might suggest that the theropod speciation pattern would be a classic “early burst” or adaptive radiation-type speciation ( 20 ) with long protracted branches ( 8 ), corresponding to a speciation slowdown model. Although our results do show an initial burst of speciation events and a gradual and prolonged slowdown, consistent with an early burst model, the fact that extinction rate surpasses speciation rate highlights a more complex process in theropods ( Figs. 2B and 3 ).

The most prominent downturn is seen in the sauropodomorphs, where speciation increases rapidly through the Triassic and Early Jurassic (an average of 0.137 speciation events for every 1 My) until ∼195 Mya; then, speciation rate starts to slow down, and extinction rate surpasses speciation rate at ∼114 Mya ( Figs. 2B and 3 ). Early sauropodomorph lineages are numerous but not long-lasting, and taxa that originated earlier in geological time are successively replaced by younger ones. The near extinction of the diplodocoids at the end of the Jurassic 145 Mya did not affect high speciation rates ( Fig. 3 ), and sauropodomorphs only began their decline ∼30 My into the Early Cretaceous ( Fig. 3 ). The subsequent originations of titanosaurian taxa were not nearly enough to compensate for the continuous loss of sauropods throughout the remainder of the Cretaceous.

Net speciation per 1 My through time in Mesozoic dinosaurs. Net speciation per 1 My can be calculated from model predictions ( Fig. 2B ) as differences between intervals (here, per 1 My). (A) Each branch of a dinosaurian phylogeny was assigned a net speciation per 1-My value based on its temporal location and group membership and plotted on a color gradient. Earlier branches have higher net speciation per 1 My (orange), whereas later branches have lower net speciation per 1 My (dark gray), except in Hadrosauriformes and Ceratopsidae, in which net speciation per 1 My increases with time. (B) The three major dinosaur groups [Sauropodomorpha (blue), Theropods (red), and nonhadrosauriform, nonceratopsid Ornithischia (green)] show an early onset of speciation slowdown until the middle of the Early Cretaceous, when speciation rates are exceeded by extinction rate (net speciation per 1 My falls below zero; dashed horizontal line). Values above zero indicate increases in species counts, whereas those below zero indicate decreases in species counts. (B, Inset) Hadrosauriforms (light green) show a slow increase in net speciation per 1 My through time, whereas ceratopsians (light blue) show a highly variable but on average, rapid increase toward the end of the Cretaceous. Posterior predictions (transparent lines) show the uncertainties in the model. Mean posterior values are shown as bold lines. Vertical lines indicate major stratigraphic boundaries (with their ages in millions of years ago) like those in Fig. 2 . Silhouettes courtesy of PhyloPic.org (CC BY 3.0) and Jack Mayer Wood (Parasaurolophus), Mathew Wedel (Brachiosaurus), Andrew A. Farke (Stegosaurus and Centrosaurus), and Martin Kevil (Carcharodontosaurus). E Cret, Early Cretaceous; L Cret, Late Cretaceous; L Triassic, Late Triassic.

Because nonhomogeneity in evolutionary rates is widespread and common in nature ( 17 ⇓ – 19 ) and dinosaurs are diverse—from the bipedal, carnivorous theropods to the quadrupedal, megaherbivorous sauropods—we might expect to find different speciation dynamics in the different dinosaurian subclades. When model parameters were estimated separately for each of the three main subclades (Ornithischia, Sauropodomorpha, and Theropoda), the same general pattern as in the total Dinosauria model was recovered but with extinction rates exceeding speciation rates earlier at 48–53 My before the K-Pg boundary (ΔDIC > 12) ( Fig. 2B and Table S1 ). Ornithischia here refers to nonhadrosauriform, nonceratopsid ornithischians, because the two Cretaceous subclades, Hadrosauriformes and Ceratopsidae, show speciation patterns distinct from other ornithischians; Lloyd et al. ( 8 ) also identified significant diversification shifts at the base of comparable clades [i.e., Euhadrosauria (here Hadrosauriformes) ( SI Text ) and Ceratopsidae]. In line with this finding, these two subclades show no signs of speciation slowdowns or downturns (ΔDIC between linear and quadratic models >5 in favor of the linear model) ( Figs. 2B , Inset and 3 and Table S1 ). Thus, the difference in the timing of the switch from slowdown to downturn in the Dinosauria model and for the three major clades is caused by the nonhomogeneity in speciation processes across dinosaurian groups. However, these two subclades combined only represent 14% of dinosaur species; over time, dinosaurs overwhelmingly experienced a reduction in their capacity to replace extinct species with new ones—net speciation per 1 My at the time that dinosaurs went extinct (66 Mya) was significantly below zero (speciation rate is less than extinction rate) ( Fig. 3B ) in the three major clades ( Table S12 )—and Hadrosauriformes and Ceratopsidae are the exceptions.

Model predictions of speciation through time in Mesozoic dinosaurs. (A) Compared with the linear model (orange), the quadratic model displaying a speciation slowdown (dark gray) substantially improves model fit (ΔDIC > 4). (B) This pattern holds true in the three major clades [Ornithischia (green), Sauropodomorpha (blue), and Theropoda (red)] and further improves model fit. (B, Inset) Model fit significantly improves when separate model parameters are estimated for the ornithischian subclades Hadrosauriformes (light green) and Ceratopsidae (light blue) from other ornithischians, but the slowdown and downturn are not observed for the two Cretaceous ornithischian subclades. Posterior predictions (transparent lines) show the uncertainties in the model. Mean posterior values are shown as bold lines. Vertical lines indicate major stratigraphic boundaries (with their ages in millions of years ago). Silhouettes courtesy of PhyloPic.org (CC BY 3.0) and Jack Mayer Wood (Parasaurolophus), Mathew Wedel (Brachiosaurus), Andrew A. Farke (Stegosaurus and Centrosaurus), and Martin Kevil (Carcharodontosaurus). E Cret, Early Cretaceous; L Cret, Late Cretaceous; L Triassic, Late Triassic.

Using a phylogenetic generalized linear mixed model (GLMM) in a Bayesian framework ( 15 ) and three recent large comprehensive dinosaur phylogenies comprising 420 ( 8 ) and 614 taxa [two trees ( 16 )], respectively, we found that the data are significantly better explained by a model, in which extinction rate exceeds speciation rate from ∼24 My before the K-Pg boundary, than the simpler alternative model [difference in deviance information criterion (ΔDIC) between linear and quadratic models >11] ( Fig. 2A and Table S1 ). Our findings are qualitatively identical across all three trees, and we report on results from one of the 614-taxon trees ( 16 ).

As an indirect measure of the influence of geography on speciation dynamics, such as segregation by geographic barriers ( 30 ), we used Mesozoic eustatic sea-level reconstructions ( 34 ) as an additional covariate in our models (mean sea-level value along each terminal branch). We also tested the ecological limit on clade diversification or the possible effects of niche saturation by adding a measure of intraclade diversity taken as the number of contemporary branches (including internal branches) for each taxon (the number of tips in time-sliced trees) ( 48 ). All data files are available in Datasets S1–S13 .

Because the fossil record has long been known to be incomplete ( 50 , 51 ), it is possible that the observed slowdown and downturn are byproducts of undersampling. This assumption would imply that there is a systematic downward bias in the phylogeny toward recent times, which would be counter to the usual expectation for poor sampling ( 50 , 51 ). Here, to test the effect of such biases, we fitted additional models with appropriate covariates, including stage-level formation counts (because formation count is widely reported to be associated with sampling bias) ( 9 , 10 , 12 , 35 , 44 , 52 , 53 ), taxon-specific formation counts (the number of formations in which a taxon is found), taxon-specific collection count (the number of fossil collections in which a taxon is represented), cladewise valid taxa counts (the known underrepresentation in the phylogeny) ( 54 ), fossil quality scores (state of preservation) ( 55 ), and body size (smaller taxa are less likely to be preserved) ( 56 ).

Model fit was assessed using deviance information criterion (DIC) and inspection of model parameter significance (using the proportion of the posterior distribution of the MCMC estimate that crosses over zero multiplied by two). We determined the best fit model as the model with the lowest DIC score and a difference in DIC score compared with that of a base model (ΔDIC) that is greater than four. In the case where multiple models had nonsignificant differences in model fit (i.e., ΔDIC < 4), we inspected the significance of model parameters and selected the model with significant covariates (i.e., nonsignificant covariates were removed).

Chains were run for 10 6 iterations, with sampling at every 1,000th iteration. We fitted a GLMM with a Poisson link to appropriately account for error structure in count data–although we discuss predicted curve shapes in log space, we did not log-transform node count for model fitting ( 49 ). MCMCglmm automatically accounts for overdispersion in the count data distribution. We used default priors (μ = 0 and V = I × 10 10 , where I is an identity matrix) for the fixed effects and parameter-expanded priors (V = 1, ν = 1, αμ = 0, and αV = 25 2 ) for the phylogenetic random effects ( 15 ).

Separate intercepts, slopes, and quadratic terms were estimated for the three major dinosaurian clades (Sauropodomorpha, Theropoda, and Ornithischia; three-group model). Lloyd et al. ( 8 ) previously identified two significant diversification shifts in the Cretaceous ornithischians at the base of the clades Euhadrosauria (here Hadrosauriformes) and Ceratopsidae, and therefore, we estimated separate model coefficients (intercepts and slopes) for these groups from other ornithischians (five-group model).

We fitted GLMMs in a Bayesian framework through MCMC using the MCMCglmm R package ( 15 ). The total number of speciation events (node count) along the phylogenetic path for each taxon was modeled as the response variable, with the corresponding path length (time elapsed from root to tip) as the main effects predictor variable—this model formulation forms the null linear model ( Fig. 1 , model A). We also fitted a speciation slowdown model with the addition of a quadratic term (time 2 ) to the main effect. Incidentally, a quadratic model can also explain the opposite case, where speciation rate increases while extinction rate remains constant. We include phylogeny as a random effect to account for shared ancestry.

We used three recent large comprehensive dinosaur phylogenies comprising 420 ( 8 ) and 614 taxa [two trees ( 16 )]. Trees were scaled according to the midpoint time of each terminal stratigraphic range ( 16 ) using the “equal” scaling method ( 47 ) implemented in the paleotree R package ( 48 ). Additionally, we scaled the trees using two alternative sets of terminal dates, the first appearance dates and the last appearance dates, to assess the effects of tree scaling on model results.

SI Text

Phylogeny. We used three large-scale dinosaur trees: the two composite trees by Benson et al. (16) with 614 taxa each and the supertree by Lloyd et al. (8) with 420 taxa. Trees were scaled according to the midpoint time of each terminal stratigraphic range [data from Benson et al. (16)] using the equal scaling method by Brusatte et al. (47) implemented in the paleotree R package (48). Additionally, we scaled the trees using two alternative sets of terminal dates, the first appearance dates (FADs) and the last appearance dates (LADs), to assess the effects of tree scaling on model results: FAD-scaled trees would have shorter path lengths, making them more likely to reduce or eliminate the speciation slowdowns and decreases, whereas LAD may enhance them. We set the root age of the tree by Benson et al. (16), including the dinosauromorphs, using the age of the oldest taxon (Asilisaurus kongwe) in the dataset plus a very small amount of time (10−05) and pruned out the nondinosaurian dinosauromorphs. The root age of Dinosauria was determined by using a more inclusive tree to avoid the problem of assigning an arbitrary root branch length, which is necessary for branch scaling when using the method by Brusatte et al. (47). The root age of Dinosauria from this pruned tree (244.3, 245.08, and 238.88 Mya for the midpoint, FAD, and LAD scaling, respectively) was then used to set the root of the supertree by Lloyd et al. (8). Model and statistical outputs for all trees are given in Tables S1–S12.

GLMMs. To model speciation dynamics and compare whether dinosaur speciation conforms to a null linear model or a slowdown model, we fitted GLMMs in a Bayesian framework through MCMC using the MCMCglmm R package (15). The total number of speciation events (node count) along the phylogenetic path for each taxon was modeled as the response variable, with the corresponding path length (time elapsed from root to tip) as the main effects predictor variable—this model formulation forms the null linear model (Fig. 1A). We also fitted a speciation slowdown model with the addition of a quadratic term (time2) to the main effects. Incidentally, a quadratic model can also explain the opposite case, where speciation rate increases while extinction rate remains constant. We further fitted a model in which the main effects predictor variable was square root-transformed (√time), which models a slowdown toward an asymptote, and compared its fit with the quadratic model. If the quadratic model performed better and we observe a downturn, then we can conclude that a slowdown model with a downturn (Fig. 1C) is a better explanation of the data than a slowdown toward an asymptote (Fig. 1B). We include phylogeny as a random effect in the form of an inverse phylogenetic variance–covariance matrix to account for shared ancestry. Speciation slowdowns have been reported for a wide range of extant phylogenies (30), but previous analyses have not distinguished the difference between a slowdown toward an asymptote (Fig. 1B) and the actual demise of a clade (Fig. 1C), in which extinction rate exceeds speciation rate at some time in the evolutionary history of the group. By sampling numerous taxa over a long period (i.e., a tree with numerous extinct taxa, such as dinosaurs) and modeling the relationship between speciation events and time at extinction given an expected model of speciation through time, it is possible to make this distinction. The number of nodes at time 0 is zero, and thus, the model should theoretically be fitted through the origin. However, our model is more conservative, in that we allow intercepts to be estimated; by estimating the intercept, we account for potential uncertainties in the basal divergences (e.g., missing taxa or incorrect dating of the root). Separate intercepts, slopes, and quadratic terms were estimated for the three major dinosaurian clades (Sauropodomorpha, Theropoda, and Ornithischia; three-group model). Lloyd et al. (8) previously identified two significant diversification shifts in the Cretaceous ornithischians at the base of the clades Euhadrosauria (here Hadrosauriformes) and Ceratopsidae, and therefore, we estimated separate model coefficients (intercepts and slopes) for these groups from other ornithischians (five-group model). We used the more inclusive clade Hadrosauriformes (least inclusive clade containing Hypacrosaurus stebingeri and Iguanodon bernissartensis or Lurdusaurus arenatus depending on the tree) instead of Euhadrosauria (Hadrosauridae and Fontllonga euhadrosaur), because basal hadrosaurs (including Iguanodon, etc.) show similar patterns in the raw data (node count ∼ time plot) with derived hadrosaurids. Chains were run for 106 iterations, with sampling at every 1,000th iteration. We fitted a Poisson GLMM to appropriately account for error structure in count data—although we discuss theoretical curve shapes in log space, we did not log-transform node count for model fitting (49). MCMCglmm automatically accounts for overdispersion in the count data distribution. We used default priors (μ = 0 and V = I × 1010, where I is an identity matrix) for the fixed effects and parameter expanded priors (V = 1, ν = 1, αμ = 0, and αV = 252) for the phylogenetic random effects (15). Model fit was assessed using DIC and inspection of model parameter significance (using p-MCMC, the proportion of the posterior distribution of the MCMC estimate that crosses over zero multiplied by two). We determined the best fit model as the model with the lowest DIC score and a difference in DIC score compared with that of a base model (ΔDIC) that is greater than four. In the case where multiple models had nonsignificant differences in model fit (i.e., ΔDIC < 4), we inspected the significance of model parameters and selected the model with significant covariates (i.e., nonsignificant covariates were removed). Compared with the null linear models, quadratic models universally had lower (and thus, better) DIC scores. Allowing for separate model parameter estimates for the three major clades (three-group model) and subdivision of Ornithischia (five-group model) further improved model fit. Although the square root and quadratic models were not significantly different for the Dinosauria as a whole (both of which were significantly better than those in the linear model), the quadratic model was generally better than the square root model in both the three- and five-group models [except in the tree by Lloyd et al. (8)]. The five-group model with no quadratic terms for Hadrosauriformes and Ceratopsidae was by far the best model (Table S1). The slowdowns and downturns are persistent across the three trees, irrespective of branch-length scaling (midpoint, FAD, and LAD). The estimated intercept models consistently had significantly better model fit (lower DIC scores) than those with intercepts forced at zero (ΔDIC > 100). This result indicates that we are predicting more speciation events at the base of the tree than would be expected given the amount of time elapsed since the root. Either the fossil record of the earliest dinosaurs is incomplete (which obviously, is very likely) or our root estimate is too young (which according to a recent study dating the origin of dinosaurs to ∼236–234 Mya, a date considerably younger than the root constraint that we use here, would be very unlikely); thus, by estimating intercepts, we allow for these uncertainties in the model. However, we get qualitatively identical results when intercepts are fixed to zero as those when intercepts are estimated. Although the five-group models typically had lower DICs than the three-group models, the coefficient for the ceratopsid time variable is not always significant, depending on the tree [supertree by Lloyd et al. (8)] or branch lengths (scaled using FAD). In the case of the supertree by Lloyd et al. (8), this nonsignificance most likely occurs owing to the small sample size of Ceratopsidae (n = 19). Estimating a common slope between Hadrosauriformes and Ceratopsidae while estimating separate intercepts can mitigate this issue of small sample size (reducing model parameters), and we can also test how the model fit compares to a model in which separate slopes are estimated. If model fit is not significantly different, then it makes no difference whether separate slopes are estimated for hadrosauriforms and ceratopsids or a common slope is estimated between the two subclades. The p-MCMC for the common slope estimate is 0.0578, and model fit does not diminish substantially (ΔDIC = 1), indicating that the common slope model is not significantly worse than the separate slope model. Estimating separate model parameters for birds and nonavian theropods does not significantly improve model fit over the base three-group model (ΔDIC < cutoff of 4). None of the bird-specific parameters are significant, indicating that birds show nothing unique in speciation dynamics.

Sampling Biases. The fossil record has long been known to be largely incomplete (50, 51). It has been estimated that fewer than 1% of the species that ever existed are preserved as fossils (57, 58). It is, thus, widely viewed that sampling biases must be taken into account so that we do not draw incorrect conclusions from incomplete data. Although it has been widely observed that biodiversity increases through time because of better fossil preservation toward the present [the “Pull of the Recent” (59)], we explicitly test the idea that the downturn is caused by increasingly poor preservation through time. If the observed slowdown and downturn were byproducts of undersampling, then the implication would be that there is a systematic downward bias in the phylogeny toward recent times (although counter to the Pull of the Recent). Here, to account for such potential biases, we fitted additional models with appropriate covariates, including stage-level formation counts (owing to the fact that formation count is widely reported to be associated with sampling bias) (9, 10, 12, 35, 44, 52, 53), taxon-specific formation counts (the number of formations in which a taxon is found), taxon-specific collection count (the number of fossil collections in which a taxon is represented), cladewise valid taxa counts (the known underrepresentation in the phylogeny) (54), fossil quality scores (state of preservation) (55), and body size (smaller taxa are less likely to be preserved) (56). Sampling bias metrics were modeled as covariates. To represent the average amount of geological bias throughout the duration of each terminal branch, a weighted average formation count along the terminal branch was calculated for each taxon assuming that numbers of fossil-bearing formations at any given geological stage fairly represent geological bias (9, 10, 35, 44). Stage-level formation counts were taken from the work by Butler et al. (35), and for each terminal branch, weighted average values were calculated by weighting proportionally to the amount of time spent in each relevant time bin. For instance, if the branch in question spanned two time bins (time bins 1 and 2), but only 40% of its length is in time bin 1 and 60% is in time bin 2, then the weighted average of the formation counts would be [(0.4 × count 1) + (0.6 × count 2)]/2. We also quantified another measure of formation count at the taxon level using occurrence data from the Paleobiology Database (PBDB)—the number of formations from which each taxon is known was counted and treated as tip data. Likewise, we counted the number of occurrences of each taxon in the database. Valid taxa counts are the numbers of known and valid dinosaur taxa tabulated at the subclade level (e.g., Maniraptora) from two sources of data, the works by Benson et al. (16) and Benton (54), for the relevant trees (16, 60). This variable represents a known degree of undersampling in the phylogeny; many valid taxa have never been included in any numerical phylogenetic analysis [compare 965 and 1,709 valid taxa to 614 and 420 tips in the trees by Benson et al. (16) and Lloyd et al. (8), respectively]. Fossil quality score ranks the state of preservation by the amount of fossil materials preserved. For the datasets by Benson et al. (16), fossil quality score is taken as whether there is a bonebed for each taxon. For the dataset by Lloyd et al. (8), fossil quality score is as defined in the work by Benton et al. (55): whether dinosaur species were represented by one or more fragments, a complete skull, a complete skeleton, or several skeletons. Because it can be argued that there is a negative bias in fossilization potential with decreasing body size (i.e., smaller dinosaurian taxa are not as frequently fossilized as larger taxa), we modeled the effects of body mass on speciation dynamics. We pruned the tree by Benson et al. (16) to those taxa with body mass estimates available (n = 313) and fitted a five-group quadratic model with log body mass as an additional covariate. If sampling bias accounted for the speciation slowdown, then we would expect the inclusion of any of these bias metrics to eliminate the significance of the quadratic term (i.e., the slowdown model would no longer be significant). Of the six sampling bias metrics, only valid taxon count was significant, but its effect is small (for every valid taxon added, node count increases by ∼0.2%), and its inclusion in the model does not affect the quadratic term. Thus, the speciation slowdown is not a product of known undersampling in the phylogeny. The fact that formation count (both stage- and taxon-specific levels), PBDB occurrence count, fossil quality scores, and body mass are not significant indicates that availability and preservability generally do not affect speciation models. It has been argued repeatedly that fluctuations in the number of fossil-bearing formations is the primary cause of observed fluctuations in the number of unique taxa through time (9, 10, 12, 35, 44, 52). However, formation count is an inappropriate measure of sampling (22, 61), primarily because formations are not directly comparable (e.g., compare the Wessex Formation of the United Kingdom, which is only found in one United Kingdom county, with the Morrison Formation of the United States, which spans at least six states and 19 US counties; PBDB, retrieved June 2, 2015). More specifically, the number of unique genera per formation ranges from 1 to 78 (mean = 4), whereas the number of fossil collections (geographic sites) per formation ranges from 1 to 228 (mean = 6.5); the numbers of available fossiliferous locations and documented occurrences greatly vary from formation to formation, spanning two orders of magnitude. With respect to rock volume, the variation spans eight orders of magnitude (55). Thus, a simple count of formations is not a fair sampling proxy (i.e., all formations are treated equally when they are not). Furthermore, the directionality (what drives what) between cumulative metrics and paleodiversity is not generally determined. Through historical study time, the number of dinosaur taxa, the number of dinosaur localities, the number of dinosaur collections, and the number of dinosaur-bearing formations increase in tandem. In synoptic studies, however, additional information is required to determine why values fluctuate through geological time. For example, low values of paleodiversity, formations count, collections, and localities could suggest a deficit of fossils (poor sampling) or simply low diversity, perhaps in the time after a mass extinction event (no evidence about sampling). This redundancy of metrics has been noted (55, 61) but widely ignored (9, 12, 35, 44, 52, 53). Similarly, just because sampling metrics are often correlated with paleodiversity does not mean that the former drives the latter. Thus, it can also be argued that such sampling metrics and paleodiversity are unrelated but independently follow similar temporal trends (i.e., spurious regression). For instance, a bifurcating phylogenetic tree will increase the number of lineages exponentially (log-linearly) through time, whereas fossil/rock availability will follow the Pull of the Recent (59) effect. Both are inherent properties of the respective data, but there is no causal relationship between them. The state of preservation (fossil quality), however, could be a more reasonable sampling proxy compared with formation count because of its taxon-specific, taphonomic, and geologic nature. Fossil quality effectively scores the potential for any speciation event to be missing along the terminal branch; a taxon with poor preservation is likely to possess a sister taxon of equal preservation potential, the fossil of which has not thus been discovered. However, fossil quality as scored here is not significant as a covariate, and it does not affect the regression parameters for time or time2. Recent fossil discoveries and statistical analyses, especially of maniraptoran taxa, have led to a consensus view of a general tendency for body size reduction in lineages leading to birds (21, 22, 62). This observed negative trend in body size evolution could potentially translate to an overrepresentation of smaller dinosaurian taxa toward the K-Pg boundary. Small size, in turn, has been argued to be a negative bias for fossilization (56), and therefore, it allegedly follows that undersampling increases through the history of Dinosauria. However, our results do not support this view; log body mass is not a significant covariate, and more importantly, it does not diminish the quadratic effects in any of the major dinosaurian clades exhibiting slowdowns/downturns (Tables S1–S10). Furthermore, size reduction through time has only been observed in a specific clade of theropods and is not present at all in ornithischians or sauropodomorphs—size reduction is not a universal phenomenon across the different dinosaur clades.

Extrinsic Controls on Speciation Dynamics. Speciation dynamics (in particular, the slowdown) have been much discussed in the context of biogeography (reviewed in ref. 30). For instance, allopatric, parapatric, and peripatric speciation are popular geographic mechanisms for speciation. We can indirectly assess the influence of geography, such as segregation by geographic barriers (30), using Mesozoic eustatic sea-level reconstructions (34) as an additional covariate in our models. We calculated the mean sea-level value along each terminal branch instead of the terminal value at the point in time of the tips, because we are interested in the average effect of sea level on speciation for the whole duration of each taxon. Furthermore, the eustatic sea-level estimates by Haq et al. (34) are highly variable at fine timescales, and a point estimate value for a specific point in time would be subject to a wide margin of error given that fossil taxa can seldom be precisely and accurately dated and also, that the sea-level curves themselves are generalized temporally. Mean sea level was significant, with a positive regression coefficient; increases in sea level are associated with increases in speciation. Rises in sea levels are associated with landward advances of coastlines and continental flooding, the biogeographical effects of which being shrinking and fragmentation of species ranges. Such separations by geographic barriers lead to reproductive isolation and subsequent speciation; the higher the sea level, the more the barriers. Horner et al. (29) observed that the emergence of transitional morphotypes coincides with marine transgressions in Late Cretaceous rocks of western North America, consistent with our finding that increasing sea levels induce speciation. However, the effect of sea level on speciation is small; for every 1-meter increase in sea level, node count increases by 0.2–0.25%. Sea level also does not diminish the time quadratic term, and therefore, it does not account for the speciation slowdown. The fact that sea level significantly affects speciation conflicts with a previous analysis, in which no association between dinosaurian taxic diversity and fluctuations in sea level was found (35). This discrepancy may be partly because of sample size; Butler et al. (35) had only 26 geological stage-level time bins, whereas we used 420 and 614 tip data at ∼100–120 distinct time points for the datasets by Lloyd et al. (8) and Benson et al. (16), respectively. Statistical significance is more difficult to establish in smaller datasets. Furthermore, stage-level tabulation of unique taxon counts may be too coarse to capture subtle biodiversity responses to environmental fluctuations. Stability in taxon counts between two consecutive time bins can be achieved through various means [for instance, that the same taxa persist through time (gross diversity remaining constant) or that a certain proportion goes extinct but is replaced by newly originated taxa (net diversity remaining constant)]. Fluctuation in sea level would not have impacted the former, but it may have affected the latter (in extinction and origination). It is not possible to differentiate between these two (if not more) forms of stability from time bin to time bin in simple tabulated time series; tabulated data do not adequately capture the underlying evolutionary processes (which are inherently phylogenetic in nature) that led to the apparent number of taxa in a succession of time bins. In addition, previous studies did not take into account shared ancestry (taxic counts even within a time bin are not independent) and phylogenetically informed speciation and extinction rates [e.g., “origination” (speciation) does not occur independent of other taxa and ancestors (taxa do not magically appear)]. We also tested the ecological limit on clade diversification or the possible effects of niche saturation. Members of the same clade are more likely to compete for similar if not the same ecological niche (32, 33), and the more numerous the number of contemporary lineages, the fewer the number of available niches. As clades expand, available ecological niches fill up until a limit is reached, and speciation will be forced to slowdown. Intraclade diversity was measured as the number of contemporary branches (including internal branches) for each taxon by counting the number of tips of time-sliced trees using the paleotree R package (48). We used the same clades as in the three- and five-group model specifications. However, cladewise lineage diversity does not explain the observed speciation slowdown—the variable is not significant, and its inclusion does not diminish the time2 effect. It is possible that our coding of contemporaneous intraclade competitors is at a taxonomic level that is too high to accurately measure niche competition. For instance, body size alone in theropods spans at least three orders of magnitude (from ∼1 to >5,000 kg), but they also vary in diet and lifestyle (e.g., compare the likely diets of the contemporary and sympatric theropods Troodon and Tyrannosaurus), and it is unlikely that such taxa would be in direct competition with each other for the same ecological niche. However, closely related taxa are more likely to be in competition with each other than distantly related taxa (e.g., theropods and sauropods), and therefore, our analyses presented here provide a general picture. More detailed competition models may provide additional insight into this matter.

Effects of Covariates on Model Fits and Parameters. Each variable (sampling metrics and extrinsic factors) was separately added as a covariate to the five-group model. Whether the addition of a covariate improved the model was determined by the significance of the parameter estimates as well as DIC improvement over the base model using a cutoff ΔDIC value of four. Because covariates may have hidden effects in the presence of other covariates, we modeled a full set of covariates and a reduced set of covariates based on the significant variables in the full model (a form of backward elimination). The model with the reduced set of covariates improves DIC substantially over the base, single-covariate, and full-covariate set models (Table S1). The significant covariates of the reduced set are sea level and valid count; the model that accounts for known underrepresentation in the phylogeny (valid count) and a potential extrinsic factor (sea level) performs better than a simple five-group quadratic model. Importantly, even with these additional covariates accounted for, the quadratic term is still significant, and a strong downturn in net speciation per 1 My is present. Speciation slowdown and downturn in dinosaurian clades cannot be sufficiently explained even by the combination of these two variables.

Effects of Branch Lengths. Although differently scaled branch lengths (midpoint, FAD, or LAD) did not affect the overall slowdown and downturn patterns in speciation, they did affect significances of predictors. This effect is generally more prominent in trees scaled using FAD than those scaled using LAD. In the supertree by Lloyd et al. (8), using LAD had the effect of eliminating the model fit improvement of the five-group model over the three-group model (ΔDIC = 0). Using FAD to scale branches instead of range midpoints has the effect of eliminating the significance of ceratopsids in the trees by Benson et al. (16) (although model fit improvement remains significant; ΔDIC > 20), because branches are shortened, thereby reducing the variability in time (variance in x is too small compared with variance in y), which is especially severe for a late-appearing clade, such as Ceratopsidae. A recent analysis on the geochronology of the Argentinian Chañares Formation, in which early dinosauromorphs are found, resulted in an age 5–10 My younger than previously thought (from ∼245–242 Mya to ∼236–234 Mya) (63), pushing forward the timing for the origin of dinosaurs and the root age of the dinosaur phylogeny. Given this new information, it is likely that the initial burst of speciation was even higher than we have recovered in our models here and that the downturn may have been even more pronounced.

Assessing the Rate of False Detections. To assess whether speciation slowdowns (and downturns) can be falsely detected by chance in trees that were grown under constant rates of speciation, we simulated 100 constant-rate birth–death trees using the R packages phytools (64) and TreeSim (65). We then fitted linear and quadratic models to each of the 100 trees using MCMCglmm (15). We compared model fit using ΔDIC as described above. The expectation is that, in only 5% of the simulated trees, the quadratic model is a better fit than the linear model. We repeated this procedure for trees simulated under two different conditions: (i) stopping the simulation when the number of extant taxa reaches N and (ii) stopping the simulation at time t. We used N = 24 [the theoretical number of extant tips N at time t given as N(t) = 2e(λ − μ)t (66) when λ = 1, μ = 0.5, and t = 5 (67)] for each of the stopping conditions. For condition 1, we first simulated 500 trees with a sufficiently large N (N = 500) using the sim.bd.taxa R function (65) and sampled from this pool of trees using the General Sampling Approach (68) with the sim.gsa.taxa R function (65) so that N = 24. All extinct tips were kept in the simulated trees. From this sample of trees, we further randomly sampled 100 trees as our final set of trees. For condition 2, we simulated 1,000 trees to t = 5 using the pbtree R function (64) resulting in a large sample of trees, which includes completely extinct trees (i.e., those that did not make it to time point t, the vast majority of which had very low total tip count n; all of the tips in a tree including extinct taxa). Therefore, we sampled 100 trees that fulfilled the criterion n ≥ 30 to gain a tree of a sufficient size for GLMM fitting. Under both conditions, in only 4 of 100 simulated trees did we find that the quadratic model had a significantly better fit compared with the linear model. These results indicate that type II errors are within expectation at a 5% level. Significant speciation downturns do not occur by chance, and our predictions in the dinosaur trees are a genuine signal of a long-term demise in which extinction rate exceeded speciation rate tens of millions of years before the K-Pg boundary.