Article | Open | Published:

# Closing the gap between palaeontological and neontological speciation and extinction rate estimates

Nature Communicationsvolume 9, Article number: 5237 (2018) | Download Citation

## Abstract

Measuring the pace at which speciation and extinction occur is fundamental to understanding the origin and evolution of biodiversity. Both the fossil record and molecular phylogenies of living species can provide independent estimates of speciation and extinction rates, but often produce strikingly divergent results. Despite its implications, the theoretical reasons for this discrepancy remain unknown. Here, we reveal a conceptual and methodological basis able to reconcile palaeontological and molecular evidence: discrepancies are driven by different implicit assumptions about the processes of speciation and species evolution in palaeontological and neontological analyses. We present the “birth-death chronospecies” model that clarifies the definition of speciation and extinction processes allowing for a coherent joint analysis of fossil and phylogenetic data. Using simulations and empirical analyses we demonstrate not only that this model explains much of the apparent incongruence between fossils and phylogenies, but that differences in rate estimates are actually informative about the prevalence of different speciation modes.

## Introduction

Understanding the drivers of biodiversity is a major, long-standing research focus in evolutionary biology1,2,3,4,5. Since changes in taxonomic diversity reflect the combined effects of speciation and extinction processes, intense effort has gone into developing methods to quantify species diversification rates. Two primary sources of data are used to estimate speciation and extinction rates: (i) phylogenetic trees of extant taxa6,7,8,9,10 and (ii) fossil occurrence data3,11,12,13. Although extant and fossil species are samples of the same underlying diversification process, a large discrepancy has been widely documented between empirical estimates of diversification rates inferred from phylogenetic versus fossil data13,14,15,16, with few exceptions, e.g. ref. 17. In particular, extinction rates estimated from empirical phylogenies are often much lower than expected given observed extinct diversity18. For example, phylogenetic estimates of diversification rates among cetaceans suggest speciation has exceeded extinction over the past 12 Myr, implying diversity has increased towards the recent19. In contrast, analyses of the cetacean fossil record indicate extinction has exceeded speciation over this same interval, and that the diversity of cetaceans was in fact once much higher than it is today15,20,21.

Discrepancies remain despite the development of increasingly realistic approaches to estimating diversification rates9,22,23,24, and efforts to explore the limitations of both stratigraphic and phylogenetic data13,18,20,25,26. Since overwhelming palaeontological evidence suggests that extinction rates close to zero are extremely improbable, some authors have questioned whether it is possible to estimate diversification rates reliably from phylogenetic data alone15,25. This situation is surprising since (i) methods used to estimate rates from fossils and phylogenies are based on the same underlying mathematical birth−death process theory27,28, and (ii) both phylogenetic and palaeontological approaches perform well under simulation conditions8,13,26. Incongruences have been attributed to biases in the data29, lack of statistical power13,30,31 and violation of underlying model assumptions20,25. However, the magnitude of the observed discrepancies remains unexplained.

Previous studies have shown that different approaches to modelling speciation (or origination) can lead to divergent estimates of macroevolutionary parameters32,33,34,35,36,37. Three speciation modes (budding, bifurcation, and anagenesis) have been described in the palaeobiological literature that can leave a signature in the fossil record, without necessarily impacting the underlying phylogenetic tree38,39 (Fig. 1). These speciation modes are assumed to result from different biological processes, and may be interpreted differently depending on the adopted species concept, e.g. morphospecies versus evolutionary species33,40,41. For instance, cladogenesis via bifurcation has been presented as the expected outcome of vicariance or allopatric speciation, whereas budding divergence has been interpreted as the result of peripatric speciation38.

Here, we examine whether these alternative speciation modes, regardless of the species-generating mechanism, are responsible for driving incongruences between palaeontological and phylogenetic estimates of diversification. We unify budding, bifurcation, anagenesis, and extinction in a single ‘birth−death chronospecies’ (BDC) process, and explore the impact of this framework through theoretical considerations, extensive simulations and empirical analyses. We provide simple mathematical formulae linking the diversification rates inferred from fossils and phylogenies, taking into account alternative speciation modes, and demonstrate that different speciation modes have a large but predictable impact on incongruent diversification rate estimates. When we apply a standard birth−death model, only three out of nine empirical clades show consistent results between the fossil record and the phylogenies of extant species. However, eight out of nine empirical data sets support the BDC model, meaning fossil and phylogenetic rate estimates converge to the same diversification parameters once we account for different speciation modes. The single dataset that did not conform to the BDC model is characterised by challenging fossil taxonomy and low sampling of extant taxa. Together, our results demonstrate that there are fundamental differences in the meaning of speciation and extinction parameters inferred from fossils and phylogenies, which alone can explain much of the rate inconsistencies observed in empirical analyses. Finally, we show that the analysis of fossils and phylogenies under the BDC model can be informative about the relative importance of different modes of speciation. The BDC process unifies the definitions of speciation and extinction parameters, reflecting the differential effects of morphological evolution with and without branching events. Thus, the BDC model paves the way for a better integration of coherent models in palaeontology and phylogenetics.

## Results and Discussion

### The BDC model

We define ‘species’ as an identifiable taxonomic unit (a lineage) that can persist through time, give rise to other species, and become extinct. Under this definition, the fossil record includes extinct and extant species that can be identified on the basis of morphological traits or other information (e.g. location, age, etc.), whereas phylogenetic data typically include extant representatives only, identified on the basis of phenotypic and/or genetic data. Following the terminology developed in previous work35,38,39,42, we describe three modes that may give rise to the origination of a new species (Fig. 1), which together result in four distinct diversification processes: (1) Cladogenesis via budding: a speciation event that gives rise to one new species. The ancestral species persists and no extinction occurs; (2) Cladogenesis via bifurcation: a speciation event that gives rise to two new species, replacing the ancestral species, which becomes extinct; (3) Anagenetic speciation: evolutionary changes along a lineage that result in the origination of one new species and the replacement (extinction) of the ancestral species; (4) Extinction without replacement: a species becomes extinct without leaving any descendants. These distinct modes of speciation and extinction may reflect different evolutionary processes38,43, in addition to variation in taxonomic practices and species definitions33,44. Although substantial debate remains regarding the interpretation and characterisation of the speciation process, and the extent to which some extant and extinct species may even be considered real44,45, we consider the definition of morphospecies in the fossil record to be something of biological significance. In particular, we follow palaeontological practice by relying on an assignment of fossils to morphospecies, which in turn form species in our BDC model.

In a standard birth−death model, branching (cladogenetic) events occur at rate λ > 0, each giving rise to one additional species, and the termination of a branch (extinction) occurs with rate μ ≥ 027. We extend this model to incorporate the possibility of alternative modes of speciation, which requires two additional parameters35,46. At each branching event, bifurcating speciation occurs with probability β [0,1], while budding speciation occurs with probability 1 − β, and anagenetic speciation occurs along each branch with rate λa ≥ 0. We call this process with four parameters the ‘birth−death chronospecies’ process.

Phylogenetic estimates of diversification rates are typically obtained from dated trees of extant taxa using the reconstructed birth−death process6,7,10. In a phylogenetic tree, each tip is considered to be a different species, and different co-existing lineages refer to different species. However, no information about species assignment to lineages through time is available from a phylogeny, i.e. typically, we do not know if a given species prior to a branching event is identical to one of the species following the speciation event (cladogenesis via budding; Fig. 1). Existing phylogenetic approaches estimate the rate of branching and the rate of extinction, i.e. the parameters λ and μ. Under these models, extinction is assumed to be a lineage termination, and may not be associated with species replacement. Thus, phylogenetic approaches implicitly assume that all speciation events occur through cladogenesis via budding and neglect the extinction of lineages in the case of bifurcating cladogenesis and anagenetic speciation.

Palaeontological estimates of diversification rates, on the other hand, are obtained from observed stratigraphic ranges (the interval between the first and last appearances of a taxon in the fossil record) or estimated species ranges (the interval between species origination and extinction)13,47. Speciation and extinction rates are based on the frequencies of range origination and termination events recorded through time, and are therefore a function of the number of (morpho)species identified in the fossil record, regardless of the mechanisms that generated them. In other words, speciation rates inferred from fossil data, hereafter referred to as λ*, necessarily quantify all events leading to the generation of new morphospecies (with or without branching) and extinction rates (hereafter μ*) reflect all instances leading to the disappearance of a morphospecies (with or without replacement). These parameters can be expressed as the combined contribution of speciation through budding, bifurcation or anagenesis, and extinction without replacement:

$$\lambda ^ \ast = \lambda (1 - \beta ) + 2\lambda \beta + \lambda_{\mathrm a}$$
(1)

and

$$\mu ^ \ast = \lambda \beta + \lambda _{\mathrm a} + \mu .$$
(2)

In order to derive the expression for λ*, we note that (i) one new species may arise via budding cladogenesis (rate λ(1 − β)), (ii) two new species may arise via bifurcating cladogenesis (rate 2βλ where the 2 acknowledges that two new species arise), and (iii) one new species may arise via anagenesis (rate λa). Similarly, μ* is the sum of (i) extinction of a species due to bifurcating cladogenesis (rate λβ), (ii) extinction of a species due to anagenesis (rate λa), and (iii) extinction of a species without replacement (rate μ). In summary, with fossil data sets typically including only temporal ranges of morphospecies, we can estimate overall speciation and extinction rates but, in the absence of a robust phylogenetic hypothesis or near-complete sampling, we cannot identify different speciation modes38,48,49.

To illustrate the impact of different speciation and extinction processes on phylogenetic and palaeontological estimates of diversification rates, we performed simulations based on the BDC process under a broad range of parameters (see Methods). The initial simulations represent an ideal scenario where all species (extinct and extant) are sampled and correctly identified, and all origination, extinction and branching times are known without error. Then, we estimate λ, μ based on the phylogenies pruned of the extinct tips and λ*, μ* based on fossil data by maximum likelihood. These simulations demonstrate that we can correctly re-estimate λ, μ based on phylogenies, and λ*, μ* based on fossil data (Fig. 2).

### Mathematical exploration of the BDC model

For given λ* and μ* from fossils, and λ and μ from the phylogeny, we now explore whether it is possible to obtain the unknown parameters λa and β. We can write Eqs. (1) and (2) as

$$\left( {\begin{array}{*{20}{c}} {\lambda ^ \ast } \\ {\mu ^ \ast } \end{array}} \right) = \left( {\begin{array}{*{20}{c}} \lambda \\ \mu \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 1 & \lambda \\ 1 & \lambda \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\lambda _{\mathrm a}} \\ \beta \end{array}} \right).$$
(3)

Thus, we have a linear system with two unknowns (β and λa) and two equations. However, since the matrix $$\left( {\begin{array}{*{20}{c}} 1 & \lambda \\ 1 & \lambda \end{array}} \right)$$ does not have full rank, either no (λa, β) or infinitely many (λa, β) solve Eq. (3). In particular, we have three scenarios (fulfilling our parameter constraints λ > 0, μ ≥ 0, λ* > 0, μ* ≥ 0, λa ≥ 0, β [0,1]):

1. A.

Exactly one solution, namely β = 0, λa = 0, exists for

$$\lambda ^ \ast - \lambda = \mu ^ \ast - \mu \,{\mathrm{and}}\,\lambda ^ \ast - \lambda = 0,$$

which is equivalent to

$$\lambda ^ \ast = \lambda ,\,\mu ^ \ast = \mu .$$
2. B.

Infinitely many solutions exist for

$$\lambda ^ \ast - \lambda = \mu ^ \ast - \mu \,{\mathrm{and}}\,\lambda ^ \ast - \lambda > 0,$$

namely, every $$\beta \in \left[ {0,{\mathrm{min}}\left\{ {1,\frac{{\lambda ^ \ast }}{\lambda } - 1} \right\}} \right]$$ with λa = λ* − λ − λβ forms a solution.

3. C.

No solution exists if

$$\lambda ^ \ast - \lambda \ne \mu ^ \ast - \mu \,{\mathrm{or}}\,\lambda ^ \ast - \lambda < 0,$$

with the first condition coming directly from the linear system, and the second from the additional requirements of λ > 0, 0 ≤ β ≤ 1 and λa ≥ 0.

On this basis, we can specify three models with parameters λ > 0, μ ≥ 0, λ* > 0, μ* ≥ 0. The ‘equal rates model’ is appropriate if Scenario A above best describes the data. This corresponds to the budding speciation model where the palaeontological and phylogenetic parameters are the same, i.e. λ* = λ, μ* = μ. This model is a special case of the BDC model. These constraints are equivalent to λ > 0, μ ≥ 0, λ* = λ, μ* ≥ λ* − λ + μ.

The ‘compatible rates model’ is appropriate if Scenario B above best describes the data. This is the full BDC model, where λ*, λ, μ, μ* are constrained such that λ* −λ = μ* − μ and λ* ≥ λ; thus, this model contains the equal rates model as a special case. The constraints can be re-written as λ > 0, μ ≥ 0, λ* ≥ λ, μ* = λ* − λ + μ (we note that this automatically implies that μ* ≥ 0).

The ‘incompatible rates model’ is appropriate if Scenario C above best describes the fossil and phylogenetic data. Under this model, the parameters λ, λ*, μ, μ* are allowed to take any value in the range λ > 0, μ ≥ 0, λ* > 0, μ* ≥ 0, without the constraints imposed by the BDC model (Eqs. (1) and (2)); thus differences in λ and λ*, as well as μ and μ*, cannot be explained by differences in speciation mode.

We determined which of the equal, compatible or incompatible rates models are supported by different simulated and empirical data sets, by estimating λ, μ, λ*, μ* simultaneously for a given phylogeny and the corresponding fossil data (see Methods). The integration of different modes of speciation and extinction into our BDC model has substantial—in some cases unexpected—implications for the interpretation of rates estimated from phylogenetic and fossil data. Below, we highlight the most important properties of the BDC model, which emerge from the mathematical formulae above and are supported by simulations.

The first property that results from the BDC model is that phylogenetic and palaeontological speciation and extinction rate estimates will only be equal if all speciation has occurred through budding. Under the BDC model, even in an ideal scenario with fully sampled and errorless data sets, speciation and extinction rates can only be equal across phylogenetic and stratigraphic inferences if all speciation events have occurred through budding and no speciation has occurred through bifurcation or anagenesis, i.e. β = 0 and λa = 0 (Scenario A above, which is a special case of the BDC model). Any instance of cladogenesis via bifurcation or anagenetic speciation will alter the rates and contribute to an apparent incongruence between the parameters estimated from phylogenetic and fossil data, as demonstrated by simulations (Fig. 2, Supplementary Table 1). Both bifurcation and anagenetic speciation introduce additional speciation and extinction events without increasing the number of terminal taxa in the phylogeny. This will result in higher speciation and extinction rates estimated from fossil data compared to the rates estimated from phylogenetic trees of extant taxa. Thus, the BDC model predicts fossil-based rates to be equal to or exceed phylogeny-based rates:

$$\lambda ^ \ast \ge \lambda \,{\mathrm{and}}\,\mu ^ \ast \ge \mu ,$$
(4)

which directly follows from the definition of λ* and μ* in Eqs. (1) and (2). Another outcome of the BDC formulation is that phylogenetic estimates of extinction equal to zero (often inferred from empirical trees) do not necessarily imply that no extinction has occurred. Indeed, μ = 0 may indicate that extinctions have occurred through bifurcation events or anagenetic replacement (which translates to μ* > 0 in the fossil record). The term λ* = λ + βλ + λa (Eq. (1)) illustrates that bifurcating and anagenetic speciations contribute similarly in determining the discrepancy between stratigraphic and phylogenetic rate estimates.

A second property emerging from the BDC model is that phylogenetic and palaeontological estimates of net diversification will be equal irrespective of speciation mode. Phylogenetic and stratigraphic speciation rates, under the BDC model, differ by the same amount as the extinction rates, even if the parameters λa and β are unknown:

$$\lambda ^ \ast - \lambda = \mu ^ \ast - \mu .$$
(5)

Thus, while stratigraphic speciation and extinction rate estimates are likely to exceed phylogenetic speciation and extinction rate estimates, the net diversification rates are predicted to be equal:

$$\lambda ^ \ast - \mu ^ \ast = \lambda - \mu .$$
(6)

Our analyses of simulated data show that the model has the power to detect these equalities (Supplementary Figure 1).

Finally, the BDC model reveals that phylogenetic and stratigraphic speciation and extinction rates are informative about speciation mode. We cannot directly estimate the rates of cladogenesis via bifurcation and anagenetic speciation; however, we can still make some important statements regarding the prevalence of alternative speciation modes. First, the dependency λa = λ* − (1 + β)λ can be re-written as λ* − λ = λa + λβ (top row of linear system in Eq. (3)), meaning the difference between stratigraphic and phylogenetic speciation rates is the sum of anagenetic and bifurcation speciation. Second, we can provide an interval for possible values of λa. Assuming that a clade diversifies under the compatible rates model, we estimate parameters λ > 0, μ ≥ 0, λ* ≥ λ, and obtain the fourth parameter from μ* = λ* − λ + μ. Based on Eq. (3), we can obtain λa given a known β (or vice versa) from the estimated parameters λ, μ, λ*. For a known β, we obtain λa = λ* − (1 + β)λ. Since β ranges in the interval $$\left[ {0,{\mathrm{min}}\left\{ {1,\frac{{\lambda ^ \ast }}{\lambda } - 1} \right\}} \right]$$ and λa ≥ 0, we obtain λa [max{0, λ* − 2λ}, λ* − λ]. Thus, based on the estimated λ and λ*, we can provide the upper bound λ* − λ for the rate of anagenetic evolution, λa, and observe that the interval width is at most λ. Finally, we can assess the importance of anagenetic speciation relative to cladogenesis via budding. Since λ* − 2λ = λa − λ(1 − β), anagenetic speciation exceeds cladogenesis via budding when λ* − 2λ > 0, whereas budding is more frequent compared to anagenetic speciation when λ* − 2λ < 0 (Supplementary Figure 2). This property is important because it shows that the combination of fossil and phylogenetic data within the context of the BDC model, even in the absence of direct estimates of λa and β, is informative about speciation mode.

### Performance and robustness of the BDC model

If differences between phylogenetic and stratigraphic rate estimates can be explained by the BDC model, then we expect to find support for the parameter constraints λ > 0, μ ≥ 0, λ* ≥ λ, μ* = λ* − λ + μ, and thus for the relationships described in Eqs. (4)–(6). Analyses of simulated data sets using maximum likelihood, where we approximated the likelihood assuming that fossil and phylogenetic data are independent (see Methods), showed that the compatible rates model was correctly identified as the best-fitting model in >99% of cases (reported at the 0.99 confidence level) against the alternatives under a wide range of parameter settings (Table 1, Supplementary Tables 1, 2). The compatible rates model was correctly rejected in >94% of cases in favour of the equal rates model when no speciation has occurred through bifurcation or anagenesis, i.e. β and λa = 0. Support for the equal rates model reduces to zero in favour of the full BDC model as the frequency of bifurcating and anagenetic speciation events increases, i.e. β > 0 and λa > 0. Both the equal and compatible rates models were correctly rejected in favour of the incompatible rates model when fossils and phylogenies were simulated under independent processes with different rates (Table 1).

Our likelihood ratio test is robust to random, incomplete taxon sampling. Even when up to 90% of the fossil taxa are missing (Table 1, Supplementary Table 3, Supplementary Figures 3, 4), the compatible rates model was still correctly favoured over the independent rates model in 100% of the simulations. While the size of the data set certainly has an impact on the statistical power of the test, our simulations show that model testing was accurate across a wide range of realistic sampling scenarios (Supplementary Tables 3, 14). The results are also robust when fossils (and therefore ranges) are sampled under a Poisson sampling process (Supplementary Table 4, Supplementary Figures 5, 6). However, support for the compatible rates model decreases in favour of the independent rates model when sampling of fossil species is highly non-uniform (Supplementary Tables 5, 6, Supplementary Figures 710).

The power of our test was high (82–100%) even when the data were simulated under some scenarios that included a substantial amount of rate heterogeneity through time, which is not explicitly accounted for in our model (Supplementary Tables 7, 8, Supplementary Figures 11, 12). This included simulations that incorporated a period of elevated branching rates (50% increase) and a tenfold increase in extinction rate. However, the accuracy of the likelihood ratio test decreased when diversification rates varied through time and also included a period during which extinction was much greater than speciation rate, in which case the BDC model was erroneously rejected in favour of the incompatible rates model for 61–86% of replicates (Table 1, Supplementary Table 8). Importantly, these cases resulted in erroneous rejection of the BDC model in favour of the independent rates model. This suggests that our approach to model testing is conservative, in that model violations, which are likely to occur in empirical data sets, will tend to artificially decrease support for the BDC model, rather than increase it. We relaxed the assumption of constant rates by implementing a Bayesian skyline version of the BDC model in which phylogenetic and fossil-based rates of speciation and extinction can vary over time (see Methods). We demonstrate using empirical data that accounting for rate variation can improve support for the BDC.

Finally, we assessed the impact of cryptic speciation (here indicating a speciation event not accompanied by recognisable phenotypic change and therefore unobserved in the fossil record35,50) on the support for our model. As the proportion of cryptic speciation increases, support for the compatible rates model decreased most frequently in favour of the equal rates model (Table 1, Supplementary Table 9, Supplementary Figures 13, 14). This is expected since undetected speciation events in the fossil record will remove the signal of additional cladogenetic and anagenetic speciation events that occur under the BDC model. In some cases, phylogenetic rate estimates may exceed stratigraphic estimates, which is in conflict with Eqs. (4)–(6), in which case the incompatible rates model was preferred.

Taken together, the results indicate that our modelling approach correctly identifies the impact of different speciation and extinction modes on diversification rates estimated using fossils and phylogenetic data. In particular, even when the discrepancies between data sets are large (e.g. λ = 0.2, μ = 0.16 versus λ* = 0.5, μ* = 0.46), our results show that the test is able to identify support for the BDC process. However, the results also highlight several factors that can reduce the power of the test and are important to consider when applying the test to empirical data sets.

### Empirical support for the BDC model

To establish empirical support for the BDC model, we analysed fossil occurrence data and dated phylogenetic trees of nine plant and animal clades, within Bayesian and maximum likelihood frameworks (see Methods). The clades are heterogeneous in terms of temporal range, size, taxon sampling, as well as their evolutionary history and ecology (Supplementary Table 14).

Phylogenetic and stratigraphic rate estimates differ substantially in most clades (Fig. 3). These differences are statistically significant in six out of nine clades, reinforcing the observation that inconsistencies between the two data types are ubiquitous among empirical data sets (Fig. 3, Supplementary Tables 1012). Of the six data sets displaying rate discrepancies, we found that four clades conform to the expectations of the compatible rates model. Thus, although stratigraphic and phylogenetic estimates of speciation and extinction rates in these clades showed large discrepancies (e.g. 3- to 18-fold rate differences in Feliformia, Fig. 3b), these can be attributed to the occurrence of bifurcating and anagenetic speciation events, without the need to invoke any potential biases in the data.

Although the exact extent of different speciation modes cannot be inferred based on our analyses, we can use the properties of the BDC model to assess the prevalence of different speciation modes (Supplementary Table 11). For the three clades that show support for the equal rates model (Ursidae, Sphenisciformes, Canidae), we can conclude that budding was the prevalent mode of speciation, and that neither anagenesis nor bifurcation have contributed substantially to species diversification within these groups. Among the four clades that show support for the BDC model with compatible rates, our estimates indicate that anagenetic speciation was as important as budding speciation in Feliformia and Cervidae (λ* − 2λ ≈ 0) (Fig. 4). In contrast, anagenetic speciation likely exceeded budding speciation in Bovidae and Cetacea (λ* − 2λ > 0). Finally, we can obtain posterior estimates of the sum of bifurcation and anagenetic rates of speciation (λa + λβ = λ* − λ), as shown in Supplementary Figure 16.

The BDC model was rejected in two data sets (ferns and corals), indicating that rate differences may not be entirely explained by different speciation modes. These data sets are characterised by much longer evolutionary histories than the other clades (i.e. hundreds of millions of years) and exhibited the greatest amount of temporal heterogeneity in both speciation and extinction rates (Supplementary Table 13, Supplementary Figure 15), which can result in spurious rejection of a constant rate BDC model (Table 1, Supplementary Table 7). We therefore re-analysed these data sets using a skyline implementation of the model, which allows for rate variation across different intervals (see Methods). Relaxing the assumption of constant rates resulted in strong support for the BDC model among ferns within each of seven time slices used in the analysis (Supplementary Table 12, Supplementary Figure 17). Both fossil and phylogenetic data supported significant rate heterogeneity through time (Fig. 5). As expected for a genus level data set, we found evidence that budding significantly exceeds anagenetic origination throughout most of the diversification history of the group (Fig. 5a).

Thus, among the nine empirical data sets tested here, only corals show evidence for a discrepancy between phylogenetic and palaeontological estimates of diversification rates, which cannot be reconciled under the BDC model (Supplementary Table 12). This likely reflects difficulties in the identification of taxonomic lineages in corals using morphological data (the only option for the ancient fossil record) in comparison with genetic and genomic data51. The phylogenetic data also appeared to have limited signal regarding the early diversification history of the clade, as reflected by the large uncertainty in rate estimates (Supplementary Figure 18). Taxonomic incongruences between modern and fossil data are likely to introduce conflicts favouring the incompatible rates model, as we have demonstrated in simulations incorporating cryptic speciation.

### Reconciling palaeontological and phylogenetic evidence

There has been an intense, recent effort to integrate data from the fields of palaeontology and phylogenetics18,52,53,54,55. This includes not only advances in methods used to estimate diversification rates13,46,56, but also models for estimating divergence times, phylogenetic relationships54,57, and modes of phenotypic evolution53,58. Birth−death processes are fundamental to almost all methodological developments in this area of research, and consequently the definition of the parameters that underlie these models is extremely important36,38.

Speciation and extinction rates inferred from fossil and phylogenetic data are usually given equivalent definitions, that is, the expected number of speciation or extinction events per lineage per unit of time5,27,28,59. We show that interpreting these quantities as equivalent parameters requires making the assumption that all species have been generated through a budding process. In fact, in a phylogenetic framework, λ should be defined as a rate of branching, while μ is the rate of lineage extinction without replacement. Conversely, in fossil-based estimates, speciation and extinction rates quantify the pace at which (morpho)species originate and go extinct, regardless of the speciation mode or if extinction terminates a lineage or not.

The results of our empirical analyses suggest that, while our model may not be sufficient to fully capture the differences between fossil and phylogenetic evidence, the impact of alternative modes of speciation may play a previously overlooked role. Our findings also offer an explanation for striking differences in estimates of average species longevity, which are much shorter for fossil data relative to phylogenetic data16, in which extinction with replacement is unaccounted for (see also ref. 33). Understanding the conceptual differences between parameters that have previously been referred to using equivalent terms is a crucial step towards an improved integration of different data sources. Thus, the formulation of a BDC process has implications for (i) interpreting speciation and extinction rates estimated using different data sets, (ii) defining predictions about the differences between phylogenetic and fossil-based inferences, and (iii) establishing a coherent framework to simultaneously analyse phylogenetic and fossil data.

Evidently, several factors can generate incompatibility between palaeontological and phylogenetic parameter estimates, and we demonstrated some of these here, e.g. substantial rate variation, non-uniform sampling and cryptic speciation. Furthermore, speciation and extinction rates may be age-dependent16,33,60, which may be an interesting aspect for future exploration of our model. Other processes are also of potential importance, such as speciation through hybridisation and reticulation, although their integration into evolutionary models is not straightforward44. Biased sampling of fossil and phylogenetic data, taxonomic inconsistencies, and dating errors will all contribute to disparities among data sets, and have each been defended as important aspects in the development of integrative models in palaeobiology. We have shown here that part of the incongruences observed between fossil and phylogenetic estimates can be attributed to a simple yet crucial conceptual difference in the meaning of speciation and extinction rates.

### Implications of the BDC model for understanding speciation

A wide range of perspectives exist regarding the definition and the nature of species and the speciation process. For example, the evolutionary species concept tends to consider all speciation (and therefore branching) events as budding40, while the Hennigian species concept tends to consider all speciation events as bifurcating61. Meanwhile, morphospecies that arose through anagenesis (i.e. speciation via replacement) are only accepted as true species by some authors42, on the basis that this process reflects phenotypic evolution along lineages rather than speciation per se38,39.

Here, we have considered species in terms of the fundamental biological units used to obtain estimates of macroevolutionary parameters44. The description of species in the fossil record relies ultimately on phenotypic traits and we perceive changes in the pace of turnover of novel traits, or morphospecies, to reflect something of biological significance, regardless of whether these changes are associated with branching events of either type or trait turnover without branching. The combined outcome of different processes generating morphospecies will be reflected in the taxonomy of extinct taxa and therefore in the estimates of speciation and extinction rates. As we have demonstrated, these estimates will be distinct from estimates obtained from molecular phylogenies, and the discrepancies are actually informative about the speciation process under the BDC model.

The relative role of different speciation modes in the evolution of species remains largely unknown, but is likely to vary across taxonomic groups, as well as geographic and environmental contexts62. Previous work has found evidence of budding being the most prevalent mode among marine invertebrates38,49, while evidence from plants suggests that anagenesis is the predominant mode43, and conflicting evidence has been found among foraminifera48,63. However, previous work has not taken advantage of the distinction between diversification rates estimated from fossils and phylogenies. Our joint analysis of fossil and phylogenetic data provided empirical evidence of budding speciation being the most prevalent process in some clades (Supplementary Table 11, Fig. 5), but also instances supporting a substantial contribution from anagenetic speciation (while the contribution of bifurcation remains elusive; Fig. 4). By formalising the definition of the parameters that are estimated from palaeontological versus phylogenetic data and the relationship between them, we have created a new opportunity with which to approach the topic of species evolution in the fossil record.

In this paper, we fitted the BDC model treating phylogenies and fossils as independent of one another, even though both data sets were generated by the same evolutionary process. Future methodological developments should aim to infer the four parameters of the BDC while explicitly taking into account the inter-dependencies between phylogenies and fossils46,55, therefore directly quantifying the contributions of all three speciation modes.

### The sixth law of palaeobiology

Recently, C. R. Marshall15 proposed five palaeobiological laws essential for understanding the evolution of the living world, which have been broadly under-appreciated by neontologists. These laws emphasise the role of extinction in evolution, and the article joins many others in questioning the reliability of diversification rates estimated from phylogenetic trees20,25, despite a substantial body of theoretical work to the contrary10,23,26. However, the fact remains that the fossil record documents historically high levels of species diversity that are not detected from phylogenetic trees. We have demonstrated that speciation and extinction rates inferred from palaeontological and phylogenetic data are expected to differ a priori, even in the absence of any biases, simply because they measure different quantities. We also showed that a phylogenetic extinction rate of zero does not imply that species are immortal, since it ignores extinction associated with replacement, either through cladogenesis or anagenesis. Fully reconciling the discrepancies between phylogeny and fossil-based estimates of diversification rates therefore requires a better understanding of the contribution of different speciation modes to the evolution (and description) of species in the fossil record, and how these processes relate to reconstructed phylogenies.

We propose a sixth law of palaeobiology that recognises the effects of different speciation modes on the estimation and interpretation of diversification rates obtained from palaeontological and phylogenetic data. This law is given by Eqs. (1) and (2), which explicitly define the parameters associated with the processes of cladogenesis via budding or bifurcation and anagenesis, and their relationship to the diversification rates estimated from phylogenies and stratigraphic ranges. Our model illustrates that differences between fossil and phylogenetic estimates of speciation and extinction are expected and ultimately informative about the prevalent mode of speciation. The predictions of the sixth law are supported by the numerical and empirical results presented in this study, and may explain numerous other contrasting findings between phylogenetic and fossil estimates. Understanding and explicitly modelling the differences between phylogenetic and fossil species concepts should be the basis of future attempts to integrate the two data types.

## Methods

### Simulations

To validate the expected relationship between λ, μ and λ*, μ*, we used simulations incorporating phylogenetic branching processes and multiple modes of (morpho)species evolution. First we simulated constant rate birth−death trees using the R package TreeSim64. Three sets of 100 tree replicates were simulated with variable turnover: with branching rate λ = 0.2 and branch extinction rate μ = 0.02, 0.1, or 0.16. We simulated trees conditioning on the number of extant tips, n = 200. The expected origin time for each set of trees is 32, 52 and 108 time units. Discrete chronospecies units were modelled by combining three alternative modes of speciation (cladogenesis via budding, cladogenesis via bifurcation and anagenesis) using the R package FossilSim (https://github.com/fossilsim/fossilsim). In a given tree, each branching event (or node) represents a budding event with probability 1 − β; thus speciation via bifurcation occurs with probability β, and anagenetic speciation occurs along each branch with rate λa. For each set of trees we generated species units using three alternative parameter combinations: β = 0 and λa = 0, β = 0.5 and λa = 0.04, or β = 0.7 and λa = 0.16. Note that when β = 0 and λa = 0, all new species occur via budding and under these circumstances we expect λ = λ* and μ = μ*. Extant species phylogenies were generated by pruning all extinct taxa from the simulated trees. Stratigraphic range data were generated using the times of origination and extinction of the chronospecies simulated under the BDC speciation model.

### Budding versus anagenetic speciation

To demonstrate that the properties of the birth−death chronospecies model can be used to determine the relative contribution of budding versus anagenetic speciation, we generated data sets with β fixed to 0.5 and λa set to 0.1, 0.2, 0.3 (λ = 0.2). We expect λ* − 2λ to be negative when budding speciation exceeds anagenetic speciation, and positive when anagenetic speciation exceeds budding speciation (see ‘Mathematical exploration of the BDC model’ for further details of model expectations; Supplementary Table 2, Supplementary Figure 2).

### Incomplete species sampling

Incomplete sampling affects many empirical phylogenetic trees and arguably all fossil data sets, erasing a proportion of the speciation and extinction events. To examine the impact of uniformly missing data on palaeontological diversification rate estimates and on support for the birth−death chronospecies model, we excluded each stratigraphic range from the initial simulated data sets prior to analysis with probability 0.1, 0.5 and 0.9 (Supplementary Table 3, Supplementary Figures 3, 4). To examine the impact of non-uniformly missing data, we excluded entirely extinct ranges only with probability 0.1, 0.5, 0.9 (Supplementary Table 5, Supplementary Figures 7, 8), and in a separate set of analysis we excluded extant ranges only with probability 0.1, 0.5, 0.9 (Supplementary Table 6, Supplementary Figures 9, 10). Note that in these experiments, we did not remove any species from the extant species phylogenies. However, removing species from the extant species phylogenies is expected to have a similar impact on the results, i.e. uniformly missing extant tips will increase the variance in parameter estimates but not necessarily reduce support for the compatible rates model, while non-uniformly missing tips may increase variance and reduce support for the compatible rates model.

### Fossil sampling

To assess the impact of fossil sampling on our model testing approach, we simulated fossil occurrence data under a Poisson process with sampling rate ψ = 0.5. Under this model the probability of sampling a given range will be a function of range duration, i.e. shorter ranges are less likely to be sampled. In the resulting fossil occurrence data sets, the start and end of ranges will be represented by first and last appearances, and thus will underestimate the total range duration. To obtain estimates of the true range durations (i.e. true origination and extinction times) we analysed the fossil occurrence data using the program PyRate (see section ‘Analysis of empirical data’ for further details). Estimated origination and extinction times were then used as input to obtain maximum likelihood estimates of λ* and μ* for model testing.

### Rate heterogeneity

Temporal variation in diversification rates is also common across clades but is not explicitly accounted for in our model. To assess the impact of temporal rate variation on the relationship between λ, μ and λ*, μ* we simulated three sets of trees under the diversification rate shift model implemented in TreeSim. To incorporate an episode of high diversification, we assigned λ = 0.3 and μ = 0 to the interval 10–20 time units. Otherwise, all other parameters settings used to generate trees were implemented as above (λ = 0.2, μ = 0.02, 0.1, or 0.16). This resulted in three sets of tree replicates with an expected origin time of 26, 34 and 52 time units. Second, we generated trees with an episode of high diversification followed by a large increase in extinction equal to speciation in the last interval. λ = 0.2 and μ = 0.02 prior to 25 time units, λ = 0.3 and μ = 0 during the interval 15–25 and λ = μ = 0.2 during the interval 0−15; expected origin time = 40. Third, we generated trees under the scenario in which extinction rate in the final interval was much greater than speciation (i.e. the clade is in decline). λ = 0.3 and μ = 0.02 prior to 25 time units, λ = 0.2 and μ = 0.01 during the interval 15−25 and λ = 0.2 and μ = 0.3 during the interval 0−15; expected origin time = 45. In all cases stratigraphic ranges were generated from the simulated trees as above, with different combinations of β and λa. Speciation and extinction rates were estimated from completely sampled data sets (Supplementary Table 7, Supplementary Figures 11, 12).

### Cryptic speciation

We define cryptic speciation as a speciation event that does not involve any recognisable phenotypic change, a case that is not accounted for in the birth−death chronospecies model. Cryptic lineages are likely to be common in clades with a fragmentary fossil record or when phenotypic changes have low probability of being preserved (e.g. traits associated with soft tissues). To examine the impact of cryptic speciation on palaeontological diversification rate estimates and on support for the birth−death chronospecies model, we generated three data sets incorporating cryptic species. At each speciation event (budding, bifurcating or anagenetic) we assigned a probability κ = 0.2, 0.5, or 0.9 of the event being cryptic, in which case the new species becomes undistinguishable from its parent species in the fossil range data set. Note this does not affect the data used to estimate rates from phylogenies. Speciation and extinction rates were estimated from completely sampled data sets (Supplementary Table 9, Supplementary Figures 13, 14).

### Incompatible rates

To establish that our test correctly rejects support for the BDC model when phylogenetic and palaeontological rates are generated under the independent rates model, we simulated data sets of trees and ranges with independent sets of parameters. We generated two sets of 100 trees under a constant rate birth−death model, in each replicate i sampling the parameter values from uniform distributions: $$\lambda _i^1\sim {\cal U}(0.1,1.5)$$ and $$\mu _i^1\sim {\cal U}(0,\lambda _i)$$ and $$\lambda _i^2\sim {\cal U}(0.1,1.5)$$ and $$\mu _i^2\sim {\cal U}(0,\lambda _i)$$. One set of trees was used to generate extant phylogenies while the second set was used to generate data sets of ranges. We then analysed sequential pairs of trees and ranges from each independent set using maximum likelihood and tested among the equal, compatible, and incompatible rate models. Anagenesis and bifurcating speciation were not incorporated into these simulations (i.e. λa = 0 and β = 0).

### Parameter inference

Phylogenies of extant species and temporal range data were used to calculate speciation and extinction rates using phylogenetic and palaeontological approaches, respectively. Maximum likelihood estimates of λ and μ were estimated using the birth−death model described by Stadler65 (Eq. (2)), which requires information about the origin time (t0), i.e. the stem age of the clade, and the age of internal nodes ($${\cal T} = [t_1,...,t_{n - 1}]$$) in a phylogeny of n extant species,

$$f({\cal T}|t_0,\lambda ,\mu ,\rho ) = \frac{{p_1(t_0)}}{{1 - p_0(t_0)}}\mathop {\prod}\limits_{i = 1}^{n - 1} \lambda p_1(t_i)$$
(7)

with,

$$\begin{array}{l}p_0(t) = 1 - \frac{{\rho (\lambda - \mu )}}{{\rho \lambda + (\lambda (1 - \rho ) - \mu )e^{ - (\lambda - \mu )t}}},\\ p_1(t) = \frac{{\rho (\lambda - \mu )^2e^{ - (\lambda - \mu )t}}}{{(\rho \lambda + (\lambda (1 - \rho ) - \mu )e^{ - (\lambda - \mu )t})^2}}.\end{array}$$
(8)

The parameter ρ is the probability of including an extant species into the phylogeny. As we included all extant species in all analyses of simulated data, we set ρ = 1 throughout. We constrained the diversification rate to be positive (i.e. λ > μ) in all analyses. We emphasise that this approach does not use any information about the delimitation of chronospecies.

Maximum likelihood estimates of λ* and μ* based on stratigraphic range data $$\left( {\cal R} \right)$$ were calculated using the birth−death model described in refs 13,28, based on the number of birth events (B) and death events (D), and on the sum of range durations (S),

$$P({\cal R}|\lambda ^ \ast ,\mu ^ \ast ) \propto \lambda ^{ \ast {B}}\mu ^{ \ast {D}}e^{ - (\lambda ^ \ast + \mu ^ \ast ){S}}.$$
(9)

In all cases we assumed complete sampling, either of extant or extinct species, unless otherwise specified.

The results using data simulated under complete sampling are shown in Supplementary Table 1 and Fig. 2, Supplementary Figure 14.

### Likelihood ratio test

To assess support for the birth−death chronospecies model under the above simulation conditions and to establish when support for the model is expected to break down, we implemented a likelihood ratio test, comparing three alternative models. Each model describes the distribution of phylogenies and the distribution of stratigraphic ranges with likelihood function:

$$P({\cal T}|\lambda ,\mu ),\quad \quad P({\cal R}|\lambda ^ \ast ,\mu ^ \ast ),$$
(10)

where $$P({\cal T}|\lambda ,\mu )$$ is defined by Eq. (7) and $$P({\cal R}|\lambda ^ \ast ,\mu ^ \ast )$$ by Eq. (8). We obtain the maximum likelihood estimates for λ, μ by maximising $$P({\cal T}|\lambda ,\mu ,)$$, and for λ*, μ* by maximising $$P({\cal R}|\lambda ^ \ast ,\mu ^ \ast )$$. We determine the model best describing our data by using a likelihood ratio test. We approximate the joint probability

$$P({\cal T},{\cal R}|\lambda ,\mu ,\lambda ^ \ast ,\mu ^ \ast ) \approx P({\cal T}|\lambda ,\mu ,)P({\cal R}|\lambda ^ \ast ,\mu ^ \ast ),$$
(11)

and then perform a likelihood ratio test. Thus we calculate LR = 2(log MLH1 − log MLH0), where $${\mathrm{ML}}_H = {\mathrm{max}}_{\lambda ,\mu ,\lambda ^ \ast ,\mu ^ \ast }P_H({\cal T},{\cal R}|\lambda ,\mu ,\lambda ^ \ast ,\mu ^ \ast )$$. Since the three models differ by the constraints imposed on parameters λ, μ, λ*, μ* we use the following rules to compare the statistical fit of alternative models.

In the Equal rates (H0) versus Compatible rates (H1) comparison the constraint λ* = λ is relaxed to λ* ≥ λ. Thus, LR = 2(log MLCRM − log MLERM) is a mixture of a χ2 distribution with one degree of freedom (i.e. a $$\chi _1^2$$) and a dirac delta distribution (i.e. all probability mass is at 0) because λ* − λ is restricted to zero in the null hypothesis and therefore lies on the border of possible values λ* − λ > 0 under the alternative hypothesis66. Thus, we reject equal rates at a level α if $$P_{\chi _1^2}(X \,> \, {\mathrm{LR}}) \,> \, 2\alpha$$. E.g., at a level α = 0.05, we reject the equal rates model if LR > 2.71.

In the Compatible rates (H0) versus Incompatible rates (H1) comparison the constraint λ* ≥ λ is relaxed to λ* > 0, and μ* = λ* − λ + μ to μ* ≥ 0. However, we only test the constraint μ* as we could not find the distribution of the likelihood ratio statistic taking into account both constraints. For this we use a $$\chi _1^2$$ distribution for the distribution of LR = 2(log MLIRM − log MLCRM), since in log-space (i.e. considering log(μ*) instead of μ*) the constraint translates to a fixed parameter versus a parameter taking any value in (−∞, ∞). Generally, the real type-1 error of our test should be slightly higher than α, as we only consider the change in the constraint on μ* and not on λ*. However, this should increase the power compared to a test at the level α. The test should therefore be conservative by slightly under-estimating the threshold yielding α = 0.05, thus favouring the rejection of the compatible rates model.

In the Equal rates (H0) versus Incompatible rates (H1) comparison the constraint λ* = λ is relaxed to λ* > 0, and μ* = λ* − λ + μ to μ* ≥ 0. We assume a $$\chi _2^2$$ distribution for LR = 2(log MLIRM − log MLERM), since for the parameters in log-space, our setting translates to two fixed parameters versus both parameters taking any value in (−∞,∞).

We note that Table 1 (main text) summarises the empirically determined type-1 errors for our likelihood ratio test in some simulation scenarios, revealing that our approximation in calculating the joint probability (Eq. (10)) and our rejection procedure produces the expected results.

Model parameters were estimated using maximum likelihood optimisation for combined data sets of phylogenies and stratigraphic ranges. All maximum likelihood optimisations were repeated five times using different initial values to reduce the probability of finding a local optima, and the results with the highest likelihood score were selected. We performed model testing using the likelihood ratio tests described above. In our tests, we used two thresholds for statistical significance set to 0.95 and 0.99.

### Analysis of empirical data

Using empirical data available for nine clades, we assessed whether there was significant incongruence between the diversification rates estimated from different data sets and, if so, whether any incongruences could be explained by the BDC model. We obtained empirical data sets from the following sources: (i−iii) Feliformia, Canidae, Ursidae, fossil occurrence data from ref. 67, phylogenies from ref. 68; (iv) Cetacea, phylogeny from ref. 69, fossil occurrence data retrieved from the Paleobiology Database (https://paleobiodb.org/) using the R library paleobioDB70; (v) Ferns and allies (genus level), phylogeny and fossil data from ref. 71; (vi—vii) Bovidae and Cervidae, phylogeny and fossil data from ref. 17; (viii) Scleractinia, phylogeny and fossil data from ref. 72; (ix) Sphenisciformes, phylogeny from ref. 73 and fossil occurrence data retrieved from the Paleobiology Database, as above. Fossil data comprise fossil occurrence times for each extinct and extant species with a known fossil record. Phylogenetic data consisted of dated phylogenetic trees of extant taxa.

To incorporate uncertainties associated with the fossil record, in addition to the maximum likelihood inference described above, we analysed the empirical data within a Bayesian framework, using a new implementation of the program PyRate74 developed for this study. This allowed us to analyse fossil and phylogenetic data to jointly estimate: (i) the times of origination and extinction of each fossil lineage, (ii) preservation rates through time, (iii) fossil-based speciation and extinction rates (λ*, μ*), and (iv) phylogenetic speciation and extinction rates (λ, μ). We estimated independent preservation rates within each geological epoch assuming a time-variable Poisson process and constant rate birth−death models for fossils and phylogeny (with independent parameters λ*, μ* and λ, μ). The birth−death processes were estimated using the likelihood functions described in Eq. (10), and we used PyRate’s default gamma priors on the birth−death rates (with shape and rate equal to 1.1). Because several of the empirical phylogenies did not include 100% of the known extant taxa, we corrected for missing lineages by setting the ρ parameter (see Eq. (7)) as specified in Supplementary Table 14. Since for the empirical phylogenies we do not know the age of the origin and we conditioned the process on the age of the crown, rather than the origin65. Although the stratigraphic and phylogenetic rates here are assumed to be fully independent parameters (as in the incompatible rates model), we used their joint posterior distributions sampled using Markov Chain Monte Carlo (MCMC) to assess the support for each model. We ran 20 million MCMC iterations to obtain posterior estimates of the parameters and used posterior samples of λ*, μ*, λ, and μ to verify the conditions predicted by the birth−death chronospecies model (Eqs. (4)–(6)).

The estimated times of origination and extinction and the phylogenies of extant taxa were used to test the equal, compatible and incompatible rates models under the maximum likelihood framework described above. Under the Bayesian framework, the equal rates model was selected if 0 was included in the 95 or 99% credible intervals of λ* − λ and μ* − μ. The compatible rates model was preferred if stratigraphic and phylogenetic rates were different, but 0 was included in the 95% (or 99%) credible interval of (λ* − λ) − (μ* − μ) and if P(λ* ≥ λ) > 0.05 (or 0.01). The incompatible rates model was preferred if none of the conditions above were met. We then re-analysed the data sets for which the BDC model was preferred, after constraining the parameter values sampled by the MCMC based on the assumptions of the equal or compatible rates models. We used the posterior samples of λ and λ* to make inferences about the prevalence of different speciation modes (Fig. 4, Supplementary Table 11).

We used the estimated times of origination and extinction for fossil lineages to infer the amount of rate variation from fossil data only. We used the reversible-jump MCMC algorithm75 implemented in PyRate to infer the number and temporal placement of rate shifts and to obtain the marginal rates through time13,76. We then computed the ratio between the greatest and the smallest marginal rates (independently for speciation and extinction) as a measure of the magnitude of rate variation in the data (Supplementary Figure 15).

In addition, we implemented a BDC skyline model in which speciation and extinction rates may vary across predefined time bins. This extension, only available within the Bayesian implementation, is based on the birth−death models with rate shifts described in ref. 22 for phylogenetic data and in ref. 13 for fossil data. In the analysis of the fern and coral data sets, time bins were set to 25 Myr in length starting from time 0 (the present) going back to 150 Ma, with the earliest bin extending from 150 Ma to the time of origin of the clade (>400 Ma). This partition scheme was selected to guarantee sufficient statistical power to estimate speciation and extinction rates with both phylogenetic and fossil data. As with the constant rate model, we first ran the analysis under the assumption of independent rates to assess whether the BDC model was supported. We then ran another analysis on the fern data under the BDC model to estimate compatible phylogenetic and fossil rates and the prevalence of different speciation modes.

### Empirical simulations

To demonstrate that we should expect to find support for the birth−death chronospecies (compatible rates) model given the scale of our empirical data (in terms of phylogenetic and stratigraphic range data size, age and diversification rates), we simulated data sets under the BDC model, based on the parameters obtained for each of the nine data sets (Supplementary Table 14). First, we estimated speciation and extinction rates from the empirical phylogenies $$\left( {\hat \lambda ,\hat \mu } \right)$$, and used them to simulate trees after setting the number of terminal tips to the present diversity of each clade. Simulated trees were then used to (1) simulate fossil ranges and (2) simulate phylogenies of extant taxa. We simulated stratigraphic range data using β = 0.5 (this value was chosen arbitrarily) and for each simulated tree, λa was iteratively increased from 0.1 until the number of simulated ranges was ≥ the empirical number of ranges. Prior to analysis the range data were uniformly pruned to match the number of empirical ranges. We generated phylogenies of extant taxa by pruning all extinct tips. For clades with incomplete taxon sampling (ρ < 1) we additionally removed a random set of taxa to reach the observed sampling fraction. Model testing using maximum likelihood was performed as described above. These simulations demonstrate that we should expect to find overall support for the birth−death chronospecies model given the scale and parameters of our empirical data sets (or similar data sets; Supplementary Table 14).

### Code availability

The maximum likelihood implementation of the BDC model (constant rate model only) is available in the R package fbdR (https://github.com/rachelwarnock/fbdR). The Bayesian implementation, which includes the BDC skyline model, is available in the latest version of PyRate74 (https://github.com/dsilvestro/PyRate). The code and scripts used in this study are also available at https://doi.org/10.5281/zenodo.1471499 (http://zenodo.org/record/1471499)77 accompanied by readme files explaining their use.

## Data availability

All the simulated and empirical data (fossil occurrences and phylogenetic trees) presented and analysed in this study are available online at https://doi.org/10.5281/zenodo.1471499 (http://zenodo.org/record/1471499)77.

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## References

1. 1.

Raup, D. M. Taxonomic survivorship curves and Van Valen’s law. Paleobiology 1, 82–96 (1975).

2. 2.

Sepkoski, J. J. A kinetic-model of phanerozoic taxonomic diversity; i. analysis of marine orders. Paleobiology 4, 223–251 (1978).

3. 3.

Alroy, J. et al. Phanerozoic trends in the global diversity of marine invertebrates. Science 321, 97–100 (2008).

4. 4.

Niklas, K. J., Tiffney, B. H. & Knoll, A. H. Patterns in vascular land plant diversification. Nature 303, 614–616 (1983).

5. 5.

Nee, S. Birth−death models in macroevolution. Annu. Rev. Ecol. Evol. Syst. 37, 1–17 (2006).

6. 6.

Nee, S., May, R. M. & Harvey, P. H. The reconstructed evolutionary process. Philos. Trans. R. Soc. B 344, 305–311 (1994).

7. 7.

Gernhard, T. The conditioned reconstructed process. J. Theor. Biol. 253, 769–778 (2008).

8. 8.

Stadler, T. On incomplete sampling under birth-death models and connections to the sampling-based coalescent. J. Theor. Biol. 261, 58–66 (2009).

9. 9.

Morlon, H., Parsons, T. L. & Plotkin, J. B. Reconciling molecular phylogenies with the fossil record. Proc. Natl. Acad. Sci. USA 108, 16327–16332 (2011).

10. 10.

Stadler, T. Recovering speciation and extinction dynamics based on phylogenies. J. Evol. Biol. 26, 1203–1219 (2013).

11. 11.

Sepkoski, J. J. Rates of speciation in the fossil record. Philos. Trans. R. Soc. B 353, 315–326 (1998).

12. 12.

Liow, L. H. & Nichols, J. D. Estimating Rates and Probabilities of Origination and Extinction Using Taxonomic Occurrence Data: Capture−Recapture Approaches, 81–94 (University of California Press, Berkeley, CA, 2010).

13. 13.

Silvestro, D., Schnitzler, J., Liow, L. H., Antonelli, A. & Salamin, N. Bayesian estimation of speciation and extinction from incomplete fossil occurrence data. Syst. Biol. 63, 349–367 (2014).

14. 14.

Herrera, J. P. Primate diversification inferred from phylogenies and fossils.Evolution 71, 2845–2857 (2017).

15. 15.

Marshall, C. R. Five paleobiological laws needed to understand the evolution of the living biota. Nat. Eco. Evol. 1, 0165 (2017).

16. 16.

Hagen, O., Andermann, T., Quental, B. T., Antonelli, A. & Silvestro, D. Estimating age-dependent extinction: contrasting evidence from fossil and phylogenies. Syst. Biol. 67, 458–474 (2018).

17. 17.

Cantalapiedra, J. L., Fernández, M. H., Azanza, B. & Morales, J. Congruent phylogenetic and fossil signatures of mammalian diversification dynamics driven by Tertiary abiotic change. Evolution 69, 2941–2953 (2015).

18. 18.

Quental, T. & Marshall, C. R. Diversity dynamics: molecular phylogenies need the fossil record. Trends Ecol. Evol. 25, 434–441 (2010).

19. 19.

Rabosky, D. L. Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS ONE 9, e89543 (2014).

20. 20.

Liow, L., Quental, T. & Marshall, C. When can decreasing diversification rates be detected with molecular phylogenies and the fossil record? Syst. Biol. 59, 646–659 (2010).

21. 21.

Pimiento, C. et al. The pliocene marine megafauna extinction and its impact on functional diversity. Nat. Ecol. Evol. 1, 1100–1106 (2017).

22. 22.

Stadler, T. Mammalian phylogeny reveals recent diversification rate shifts. Proc. Natl. Acad. Sci. USA 108, 6187–6192 (2011).

23. 23.

Etienne, R. S. et al. Diversity-dependence brings molecular phylogenies closer to agreement with the fossil record. Proc. R. Soc. Lond. B 279, 1300–1309 (2012).

24. 24.

Alexander, H. K., Lambert, A. & Stadler, T. Quantifying age-dependent extinction from species phylogenies. Syst. Biol. 65, 35–50 (2015).

25. 25.

Rabosky, D. L. Extinction rates should not be estimated from molecular phylogenies. Evolution 64, 1816–1824 (2010).

26. 26.

Beaulieu, J. M. & O’Meara, B. C. Extinction can be estimated from moderately sized molecular phylogenies. Evolution 69, 1036–1043 (2015).

27. 27.

Kendall, D. G. On the generalized birth-and-death process. Ann. Math. Stat. 19, 1–15 (1948).

28. 28.

Keiding, N. Maximum likelihood estimation in the birth-death process. Ann. Stat. 3, 363–372 (1975).

29. 29.

Raia, P. et al. Rapid action in the Palaeogene, the relationship between phenotypic and taxonomic diversification in Coenozoic mammals. Proc. R. Soc. B 280, 20122244 (2013).

30. 30.

Rabosky, D. L. & Lovette, I. J. Explosive evolutionary radiations: decreasing speciation or increasing extinction through time? Evolution 62, 1866–1875 (2008).

31. 31.

Bokma, F. Problems detecting density-dependent diversification on phylogenies. Proc. R. Soc. Lond. B 276, 993–994 (2009).

32. 32.

Smith, A. B. & Patterson, C. The influence of taxonomy on the perception of patterns of evolution. Evol. Biol. 23, 127–216 (1988).

33. 33.

Ezard, T. H. G., Pearson, O. N., Aze, T. & Purvis, A. The meaning of birth and death (in macroevolutionary birth–death models). Biol. Lett. 8, 139–142 (2012).

34. 34.

Ezard, T. H. G., Thomas, G. H. & Purvis, A. Inclusion of a near-complete fossil record reveals speciation-related molecular evolution. Methods Ecol. Evol. 4, 745–753 (2013).

35. 35.

Bapst, D. W. When can clades be potentially resolved with morphology? PLoS ONE 8, e62312 (2013).

36. 36.

Huang, D., Goldberg, E. E. & Roy, K. Fossils, phylogenies, and the challenge of preserving evolutionary history in the face of anthropogenic extinctions. Proc. Natl. Acad. Sci. USA 112, 4909–4914 (2015).

37. 37.

Faurby, S., Eiserhardt, W. L. & Svenning, J.-C. Strong effects of variation in taxonomic opinion on diversification analyses. Methods Ecol. Evol. 7, 4–13 (2016).

38. 38.

Wagner, P. J. & Erwin, D. H. Phylogenetic Patterns as Tests of Speciation Models 87−122 (Columbia University Press, New York, 1995). .

39. 39.

Foote, M. On the probability of ancestors in the fossil record. Paleobiology 22, 141–151 (1996).

40. 40.

Simpson, G. G. The species concept. Evolution 5, 285–298 (1951).

41. 41.

Wheeler, Q. D. & Meier, R. Species Concepts and Phylogenetic Theory: A Debate (Columbia University Press, New York, 2000).

42. 42.

Gingerich, P. D., Cracraft, J. & Eldredge, N. The stratophenetic approach to phylogeny reconstruction in vertebrate paleontology. Phylogenetic Anal. Paleontol. 41, 77 (1979).

43. 43.

Stuessy, T. F., Crawford, D. J. & Marticorena, C. Patterns of phylogeny in the endemic vascular flora of the Juan Fernandez Islands, Chile. Syst. Bot. 15, 338–346 (1990).

44. 44.

Freudenstein, J. V., Broe, M. B., Folk, R. A. & Sinn, B. T. Biodiversity and the species concept—lineages are not enough. Syst. Biol. 66, 644–656 (2017).

45. 45.

De Queiroz, K. Species concepts and species delimitation. Syst. Biol. 56, 879–886 (2007).

46. 46.

Stadler, T., Gavryushkina, A., Warnock, R. C. M., Drummond, A. J. & Heath, T. A. The fossilized birth-death model for the analysis of stratigraphic range data under different speciation modes. J. Theor. Biol. 447, 41–55 (2018).

47. 47.

Foote, M. Inferring temporal patterns of preservation, origination, and extinction from taxonomic survivorship analysis. Paleobiology 27, 602–630 (2001).

48. 48.

Aze, T. et al. A phylogeny of Cenozoic macroperforate planktonic foraminifera from fossil data. Biol. Rev. 86, 900–927 (2011).

49. 49.

Bapst, D. W. & Hopkins, M. J. Comparing cal3 and other a posteriori time-scaling approaches in a case study with the pterocephaliid trilobites. Paleobiology 43, 49–67 (2017).

50. 50.

Patzkowsky, M. E. A hierarchical branching model of evolutionary radiations. Paleobiology 21, 440–460 (1995).

51. 51.

Fukami, H. et al. Conventional taxonomy obscures deep divergence between pacific and atlantic corals. Nature 427, 832–835 (2004).

52. 52.

Fritz, S. et al. Diversity in time and space: wanted dead and alive. Trends Ecol. Evol. 9, 509–516 (2013).

53. 53.

Slater, G. & Harmon, L. J. Unifying fossils and phylogenies for comparative analyses of diversification and trait evolution. Methods Ecol. Evol. 4, 699–702 (2013).

54. 54.

Heath, T. A., Hulsenbeck, J. P. & Stadler, T. The fossilized birth-death process for coherent calibration of divergence-time estimates. Proc. Natl. Acad. Sci. USA 111, 2957–2966 (2014).

55. 55.

Pennell, M. W. & Harmon, L. J. An integrative view of phylogenetic comparative methods: connections to population genetics, community ecology, and paleobiology. Ann. N.Y. Acad. Sci. 1289, 90–105 (2013).

56. 56.

Didier, G., Fau, M. & Laurin, M. Likelihood of tree topologies with fossils and diversification rate estimation. Syst. Biol. 66, 964–987 (2017).

57. 57.

Ronquist, F. et al. A total-evidence approach to dating with fossils, applied to the early radiation of the Hymenoptera. Syst. Biol. 61, 973–999 (2012).

58. 58.

Silvestro, D. et al. Early arrival and climatically-linked geographic expansion of New World monkeys from tiny African ancestors. Syst. Biol. https://doi.org/10.1093/sysbio/syy046 (2018).

59. 59.

Foote, M. Origination and extinction components of taxonomic diversity: General problems. Paleobiology 26, 74–102 (2000).

60. 60.

Hagen, O., Hartmann, K., Steel, M. & Stadler, T. Age-dependent speciation can explain the shape of empirical phylogenies. Syst. Biol. 64, 432–440 (2015).

61. 61.

Hennig, W. Phylogenetic Systematics (University of Illinois Press, Urbana, 1966).

62. 62.

Benton, M. J. & Pearson, P. N. Speciation in the fossil record. Trends Eco Evol. 16, 405–411 (2001).

63. 63.

Strotz, L. C. & Allen, A. P. Assessing the role of cladogenesis in macroevolution by integrating fossil and molecular evidence. Proc. Natl. Acad. Sci. USA 110, 2904–2909 (2013).

64. 64.

Stadler, T. Simulating trees with a fixed number of extant species. Syst. Biol. 60, 676–684 (2011).

65. 65.

Stadler, T. How can we improve accuracy of macroevolutionary rate estimates? Syst. Biol. 62, 321–329 (2012).

66. 66.

Chernoff, H. On the distribution of the likelihood ratio. Ann. Math. Statist. 25, 573–578 (1954).

67. 67.

Pires, M. M., Silvestro, D. & Quental, T. B. Continental faunal exchange and the asymmetrical radiation of carnivores. Proc. Roy. Soc. B 282, 20151952 (2015).

68. 68.

Nyakatura, K. & Bininda-Emonds, O. R. P. Updating the evolutionary history of Carnivora (mammalia): a new species-level supertree complete with divergence time estimates. Bmc Biol. 10, 12 (2012).

69. 69.

Steeman, M. E. et al. Radiation of extant cetaceans driven by restructuring of the oceans. Syst. Biol. 58, 573–585 (2009).

70. 70.

Varela, S. et al. paleobioDB: an R package for downloading, visualizing and processing data from the Paleobiology Database. Ecography 38, 419–425 (2015).

71. 71.

Lehtonen, S. et al. Environmentally driven extinction and opportunistic origination explain fern diversification patterns. Sci. Rep. 7, 4831 (2017).

72. 72.

Simpson, C., Kiessling, W., Mewis, H., Baron-Szabo, R. C. & Müller, J. Evolutionary diversification of reef corals: a comparison of the molecular and fossil records. Evolution 65, 3274–3284 (2011).

73. 73.

Gavryushkina, A. et al. Bayesian total-evidence dating reveals the recent crown radiation of penguins. Syst. Biol. 66, 57–73 (2017).

74. 74.

Silvestro, D., Salamin, N. & Schnitzler, J. PyRate: a new program to estimate speciation and extinction rates from incomplete fossil record. Methods Ecol. Evol. 5, 1126–1131 (2014).

75. 75.

Green, P. J. Reversible jump Markov chain Monte Carlo and Bayesian model determination. Biometrika 82, 711–732 (1995).

76. 76.

Silvestro, D., Salamin, N., Antonelli, A. & Meyer, X. Improved estimation of macroevolutionary rates from fossil data using a Bayesian framework. Preprint at https://doi.org/10.1101/316992 (2018).

77. 77.

Silvestro, D., Warnock, R. C. M., Gavryushkina, A. & Stadler, T. Data and code http://zenodo.org/record/1471499 (2018).

## Acknowledgements

We thank David Bapst and Charles Marshall for providing valuable feedback on the manuscript. In addition, we thank Ziheng Yang, Andreas Steingötter and the ETH Seminar for Statistics group for providing advice on model testing and Carl Simpson for providing data. D.S. received funding from the Swedish Research Council (2015-04748) and from the Swedish Foundation for Strategic Research. R.C.M.W. was funded by the ETH Zürich Postdoctoral Fellowship and Marie Curie Actions for People COFUND programme. T.S. is supported in part by the European Research Council under the Seventh Framework Programme of the European Commission (PhyPD: grant agreement number 335529). A.G. was funded by the Bioprotection Research Centre. Part of the analyses were run at the high-performance computing centre Vital-IT of the Swiss Institute of Bioinformatics (Lausanne, Switzerland).

## Author information

### Author notes

1. These authors contributed equally: Daniele Silvestro, Rachel C. M. Warnock.

### Affiliations

1. #### Department of Biological and Environmental Sciences, University of Gothenburg, 41319, Gothenburg, Sweden

• Daniele Silvestro
2. #### Global Gothenburg Biodiversity Centre, 41319, Gothenburg, Sweden

• Daniele Silvestro
3. #### Department of Computational Biology, University of Lausanne, Lausanne, 1015, Switzerland

• Daniele Silvestro
4. #### Swiss Institute of Bioinformatics (SIB), 1015, Lausanne, Switzerland

• Daniele Silvestro
• , Rachel C. M. Warnock
•  & Tanja Stadler
5. #### Department of Biosystems Science & Engineering, Eidgenössische Technische Hochschule Zürich, 4058, Basel, Switzerland

• Rachel C. M. Warnock
•  & Tanja Stadler
6. #### Department of Biochemistry, University of Otago, Dunedin, 9054, New Zealand

• Alexandra Gavryushkina

### Contributions

D.S., R.C.M.W., A.G., and T.S. designed the study and developed the methods. D.S. and R.C.M.W. wrote the manuscript with contributions from all authors. R.C.M.W. implemented the maximum likelihood tests and ran the simulations. D.S. implemented the Bayesian version of the BDC model and analysed the empirical data sets.

### Competing interests

The authors declare no competing interests.

### Corresponding author

Correspondence to Daniele Silvestro.
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 2, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.8777344226837158, "perplexity": 3343.9621303872823}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.3, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2019-13/segments/1552912203547.62/warc/CC-MAIN-20190325010547-20190325032547-00215.warc.gz"}