Abstract
The branching times of molecular phylogenies allow us to infer speciation and extinction dynamics even when fossils are absent. Troublingly, phylogenetic approaches usually return estimates of zero extinction, conflicting with fossil evidence. Phylogenies and fossils do agree, however, that there are often limits to diversity. Here, we present a general approach to evaluate the likelihood of a phylogeny under a model that accommodates diversitydependence and extinction. We find, by likelihood maximization, that extinction is estimated most precisely if the rate of increase in the number of lineages in the phylogeny saturates towards the present or first decreases and then increases. We demonstrate the utility and limits of our approach by applying it to the phylogenies for two cases where a fossil record exists (Cetacea and Cenozoic macroperforate planktonic foraminifera) and to three radiations lacking fossil evidence (Dendroica, Plethodon and Heliconius). We propose that the diversitydependence model with extinction be used as the standard model for macroevolutionary dynamics because of its biological realism and flexibility.
1. Introduction
The fossil record tells us that the fate of most species is to go extinct [1]. However, it tells an incomplete story, because most species leave no trace in the fossil record. Over the last two decades, molecular phylogenies of extant taxa have come to the fore as an additional source of information about diversification dynamics, and are especially valuable for clades with a poor fossil record [2,3]. Even though extinct species are necessarily absent from a molecular phylogeny of extant species, if perlineage rates of speciation and extinction have been constant through time, extinction leaves a characteristic signature in the phylogenetic branching pattern, known as the pullofthepresent [4,5]. On a lineagesthroughtime (LTT) plot—a semilogarithmic plot with time before present on the xaxis and number of thenextant lineages with extant descendants on the yaxis—the pullofthepresent is seen as a steepening of the slope of lineage accumulation towards the present. With the number of specieslevel molecular phylogenies exploding over the last decade, it has become apparent that very few show a pullofthepresent. In fact, the opposite pattern is typical: the slope of lineage accumulation decreases towards the present, which leads to estimates of zero extinction [6–8]. Therefore, molecular phylogenies conflict with the fossil record regarding the importance of extinction, leading some to query the ability of phylogenies of extant taxa to deliver useful macroevolutionary insights in isolation from fossil data [9].
Among the hypotheses put forward to explain the puzzle of not detecting extinction in molecular phylogenies are variation of diversification rates among lineages [10], a failure to recognize the youngest species [7] and relaxation of the assumption that speciation is instantaneous [11]. Another possible explanation is that speciation rates change through time [10,12], thereby erasing or concealing the signature of extinction. While the fossil record and molecular phylogenies provide conflicting estimates of extinction rate, both support the existence of a negative feedback of diversity on the diversification rates, termed diversitydependence. Palaeontological evidence for diversitydependence comes from the relatively constant diversity within numerous higher taxa over millions of years [13–18]. Molecular phylogenetic evidence lies in the abovementioned decreasing rate of branching through time for many clades, manifest as a levellingoff on a LTT plot [19–21].
If diversitydependence and nonzero extinction are the rule, then models of diversification should include both processes. However, existing models that incorporate diversitydependence assume extinction to be zero [22]. This is in part because simulations reveal that extinction erases the signature of a reduction in the speciation rate through time [23–25], leading to the suggestion that extinction rates have probably been low in clades that show a slowdown in branching rate [22,26]. An additional motivation for ignoring extinction has been expediency [27]: it is relatively simple to evaluate the likelihood of a diversitydependent speciation model with no extinction (if there are no missing extant species), because such a model can be easily formulated in terms of a purebirth model with a timedependent speciation rate for which exact likelihood formulae exist [22]. With nonzero extinction, this mathematical simplification no longer holds because historical diversification rates depend on species that may have gone extinct and are therefore not observable in the phylogeny.
In this paper, we use a hidden Markov model (HMM) approach to numerically compute the likelihood of a phylogeny under a large variety of diversitydependent birth–death models of diversification (see box 1 and electronic supplementary material for details). In contrast to previous approaches (for an exception, see [28]), our method also allows us to take into account incomplete sampling of species (which occurs for two out of the five case studies in this paper; table 2) and presence of other species that have gone extinct but affected diversification rates when the crown group began to radiate. Furthermore, it enables efficient computation of the distribution of both the number of ancestral lineages of extant species and the total historical diversity conditioned on the molecular phylogeny, at each time between the stem or crown age of the clade and the present. Both the likelihood and the distribution of the total historical diversity can incorporate fossil data. In this paper, we use fossil data to (i) identify plausible extinction rates and (ii) conduct an a posteriori test of the degree to which parameters estimated from phylogenies of extant species captures diversity through time.
Computation of the likelihood under a diversitydependent birth–death model.
We compute the likelihood of a phylogeny under a general diversitydependent birth–death model, including the model of the main text. For a phylogeny with q extant species, we denote the branching times by t_{k} so that t_{1} < t_{2} < …< t_{q}_{−1}, with present time t_{p} > t_{q} _{− 1} and crown age t_{c} = t_{p} − t_{1} (see electronic supplementary material, figure S1). Mathematically, the diversification process is described by the master equation. This is a dynamical equation for the probability P_{n}(t) of having n species at time t, B 1
The first two terms correspond to transitions increasing the probability P_{n}(t): an extinction event from n + 1 to n species (first term) or a speciation event from n − 1 to n species (second term). The last term corresponds to transitions decreasing the probability P_{n}(t): an extinction event from n to n − 1 species or a speciation event from n to n +1 species.
We modify the master equation (B1) to guarantee that the diversification process leads to the observed phylogeny. Consider a time t between the branching times t_{k} _{− 1} and t_{k} so that the phylogeny has k branches. Denote by Q_{n}(t) the probability that a realization of the diversification process is consistent with the phylogeny up to time t and has n species at time t. The dynamical equation is B 2
The difference with the master equation (1) is in the first two terms. In the first term, we exclude extinction events in k species because none of the k branches in the phylogeny should become extinct. However, the species in these k branches can speciate. If that happens, either of the two daughter species can be included in the phylogeny, giving a factor 2k instead of k. Taken together with the speciation events in the n − k − 1 other species, we get the factor n + k − 1 of the second term in (B 2).
The following algorithm computes the probability of a phylogeny with q extant species:

— Initialize Q_{n}(t) at the first branching time t_{1} with Q_{n}(t_{1}) = 1 for n = 2 and Q_{n}(t_{1}) = 0 otherwise.

— For k =2, 3, …, q − 1 do

(i) Integrate (B2) from t_{k} _{− 1} to the next branching time t_{k}.

(ii) Multiply Q_{n}(t) by k λ_{n} at the branching event at time t_{k}.


— Integrate (B2) with k = q from t_{q} _{− 1} to the present time t_{p}.

— Extract component Q_{q}(t_{p}) from the probability vector Q_{n}(t).
This is the likelihood unconditioned on survival of the two crown lineages. We use a similar approach to compute the probability that the two lineages initiated at branching time t_{1} have descendants at present time, t_{p}. We define the likelihood as the probability of the phylogeny conditioned on the survival of the two lineages at t_{1}. Our computational approach is quite flexible and can be extended to the computation of the number of species through time (figures 2g,h and 3j–l) or the inclusion of extant species missing in the phylogeny (as is the case for Cetacea and Heliconius; table 2). We refer to the electronic supplementary material for further details.
After introducing our model, we first study under what conditions precise inference can be made. Then, to illustrate our method, we apply a model with diversitydependent speciation to two clades for which a fossil record exists and to three clades lacking a fossil record. The latter three clades are an arbitrary selection of wellsampled clades that have previously been subjected to highprofile studies on diversitydependence, with our chief selection criterion being that we wanted to have representatives of a range of different taxa. We thus aimed to provide an unbiased outlook on the prospects for estimating extinction from phylogenies.
2. Material and methods
(a) Model
The diversitydependent model that we consider is a straightforward extension of the constantrate (CR) birth–death model [29], where the speciation rate depends linearly on the number of species, n, as in a logistic population growth model: so that speciation and extinction rates equal each other at the ‘carrying capacity’ K, where λ_{K} = μ_{K}. This model is mathematically equivalent to a model with under the substitution K′ = λ_{0}K/(λ_{0} − μ). Here, K′ might be interpreted as the maximum number of niches that the species in the clade can occupy [30]. The parameter λ_{0} denotes the initial speciation rate (formally, when diversity is 0). This model (in either variant) is a generalization of the linear diversitydependence (DDL) model without extinction in Rabosky & Lovette [22], and we therefore call it the DDL+E model, and we change the name of Rabosky and Lovette's model to DDLE. We stress that our method to compute the likelihood is not limited to this model. It can be applied to other diversitydependent models, for example, the generalization of the exponential diversitydependence (DDX) model [22], or to models with diversitydependent extinction rates. DDX speciation seems to be favoured in birds [2,22], while diversitydependent extinction has received little empirical support from phylogenetic analyses [23,30] (but see [15]).
(b) Likelihood computation
When diversification rates are diversitydependent, species other than those occurring in the phylogeny—namely the extinct and nonsampled species—also contribute to historical diversification rates (because the speciation rate is a function of the number of species in existence at each point in time). Despite their contribution to diversification dynamics, extinct and nonsampled species are hidden from observation. This situation requires the use of a HMM approach, because the system being modelled is a Markov process, but with unobserved states. In box 1, we outline how such an approach allows us to (numerically) compute the likelihood of a phylogeny under a wide variety of diversitydependent birth–death models of diversification. See the electronic supplementary material for more details. Matlab code to implement the method (which is convertible to C or a standalone executable for Windows or Linux) is available from the corresponding author, and within the R package TreePar v. 2.2 [31] on CRAN (http://cran.rproject.org/web/packages/TreePar/index.html).
(c) Model fitting
We used Matlab's builtin optimization routine fminsearch to find the parameters that maximize the likelihood, choosing various initial conditions to minimize the risk of obtaining only a local maximum of the likelihood. The expected LTT plots and expected diversity through time conditional on the phylogeny were computed using the same mathematical procedures as for computing the likelihood itself. See the electronic supplementary material for details.
(d) Simulations
For a better understanding of the behaviour of the model, we plotted the expected number of lineages versus time (LTT) for various values of model parameters (λ_{0} = 0.8, K = 40, μ = 0, 0.1, 0.2 and 0.4) and crown ages (5, 10 and 15 Myr). The resulting 12 parameter combinations are sufficient to show the behaviour of the model: varying λ_{0} does not add a new dimension, because it is equivalent to rescaling time, the effect of which we already study by looking at different crown ages. In fact, the compound parameter of interest is the ratio of λ_{0} and μ, which we vary by changing μ. Varying K would mainly affect the time at which diversitydependence becomes important, so it is mostly a matter of timescale as well. To assess the performance of our approach, we performed 100 simulations of the DDL+E model for each of the 12 parameter combinations and then estimated the maximumlikelihood parameters from the branching times. If the mean or median of a parameter estimated across these 100 simulations deviates significantly from the input value, this indicates bias. If the range of estimates is wide, this indicates low precision.
(e) Phylogenetic and fossil data
The foraminifera phylogeny was constructed exclusively from palaeontological data and includes all known Cenozoic species—extant and extinct—within the macroperforate clade [32]. Reconstruction of such a phylogeny is made possible by the unparalleled fossil record of this clade: a conservative estimate of the specieslevel completeness of the record is that on average, a species is recovered from over 80 per cent of the 1 Myr bins during its existence [15]. We used Aze et al.'s [32] phylogeny of evolutionary species, rather than one of morphospecies, to minimize the influence of pseudoextinction and pseudospeciation (anagenetic change of one morphospecies into another). We applied our method to the branching times of extant species only. We then simulated diversification under the estimated model parameters and compared the simulated diversity through time with the number of evolutionary lineages present at each time (including those that later go extinct).
The maximum clade credibility phylogeny of Cetacea includes 87 of 89 extant species, was based on sequence data for six mitochondrial DNA genes and nine nuclear genes and used palaeontological age constraints on several nodes [33]. We used Quental & Marshall's [9] estimates, based on fossil genera [34], of species diversity through time. Quental & Marshall gave both a liberal and conservative estimate of species diversity. The liberal estimate was obtained by counting the numbers of genera sampled in a particular time period, assuming that all genera present in each time period coexisted. The conservative estimate included only genera that crossed specific timeboundaries. In both cases, species diversity was extrapolated after correcting for the incompleteness of the fossil record at the genus level and the present day ratio of species to genera.
The maximum clade credibility molecular phylogeny of Dendroica warblers included all 25 continental species in the genus, including four species that were formerly placed in other genera [22]. The tree was scaled, such that the root was placed at 5 Myr [35]. For Plethodon, we restricted our focus to the chronogram of the glutinosus clade comprising 30 species [36]. The Heliconius butterfly molecular phylogeny included 38 of 44 species in the clade and branch lengths were estimated using a ratesmoothing local clock approach [37].
3. Results
(a) Bias and precision: simulations
When clades are simulated over short timescales, extinction and diversitydependence make little difference to the shape of LTT plots (figure 1a; though parameter estimates differ among the models). When the models are simulated for longer, however, the LTT plots have a remarkably different shape (figure 1b,c), either saturating to a plateau when extinction is low, or showing an inverted Sshape when extinction is higher. In the latter case, extinction causes the initial tendency to saturate to eventually be counteracted by the ‘pullofthepresent’ [23,24].
We found that parameter estimates across simulated phylogenies were most biased and least precise for younger crown ages and higher extinction rates (table 1). Only when the data show a pattern similar to the curves in the rightmost panel of figure 1—either a clear pattern of saturation or a clear inverted Sshape—are the estimated parameters fairly reliable. Of our empirical datasets, only the foraminifera and Dendroica satisfy these criteria. Therefore, the maximumlikelihood parameter estimates for the other three groups should be interpreted with caution. Bias in parameter estimates (obtained by any method and not just maximum likelihood) particularly occurs for small sample sizes. The fact that the younger clades in our simulations generally have fewer species than older clades may therefore be responsible for the estimation bias we observe. The bias can here be understood by recognizing that, while a certain distribution of branching times may be a likely outcome under one set of parameters, the exact same pattern may be still more likely to arise under an altogether different combination of parameters (see electronic supplementary material for more details).
(b) Case studies: phylogenies for clades with a fossil record
To illustrate the utility of our approach, we first consider the Cenozoic macroperforate planktonic foraminifera, which hold a unique position in the study of macroevolution, owing to the existence of both an excellent fossil record and complete phylogeny for evolutionary lineages [32]. We find that our model with diversitydependent speciation plus diversityindependent extinction (DDL+E) provides an excellent fit to the inverted Sshaped LTTplot for the phylogeny of extant species (figure 2a). Extinction rates are clearly nonzero (figure 2c,e and table 2) and very closely match estimates of μ = 0.0979 from the fossil record [15]. Moreover, the number of species through time predicted on the basis of the maximumlikelihood model parameters matches the fossil evidence reasonably closely (figure 2g). In contrast, although the CR birth–death model allows for nonzero extinction, it fits poorly, although better than a diversitydependent model without extinction (table 2; [15]).
Cetacea (whales and dolphins) provide a second case study for which both an almost complete molecular phylogeny [33] and estimates of genus diversity through time, based on the fossil record [9], are available. The Cetacean LTT plot is almost linear (figure 2b), meaning that the branching pattern will contain little information to distinguish among a wide range of diversification models (figure 1a). While our method identifies a maximumlikelihood extinction rate that is nonzero (table 2), it is substantially lower than an average genus extinction rate of 0.11 calculated using fossils [9]. Although estimates of fossil pergenus extinction rates do not necessarily translate directly to species extinction rates [38,39], the genus extinction rate may give a lower bound estimate of the species extinction rate [40]. Repeating the likelihood calculations with μ fixed at 0.11, we find that the likelihood of the diversitydependent model (DDL+E) is only slightly lower (figure 2d) than when μ is free to vary and much higher than of the CR birth–death model. The clade is not fully saturated, as K = 227.6, whereas the number of extant species is 89. A recent study of Cetacean diversification based on phylogenetic information found evidence for elevated speciation rates during periods of ocean restructuring (between 35–31 and 13–4 Myr) [33]. With μ fixed at 0.11, we find the ocean restructuring (OR) model and DDL+E model to have a similar likelihood (table 2). In addition, the maximumlikelihood DDL+E parameters yielded a good agreement between the predicted historical diversity and that calculated from the fossil record of Cetacean genera [9] (figure 2h), demonstrating the benefits of combining fossil and phylogenetic data in terms of mutual illumination.
(c) Case studies: phylogenies for clades lacking a fossil record
We now confront our model with the phylogenies for Dendroica warblers in North America [22], Plethodon salamanders in North America [36] and Heliconius butterflies in the Neotropics [37]. Each phylogeny represents a radiation of at least 20 species and is wellresolved, but each group lacks a fossil record. Dendroica and Heliconius are sometimes considered to be adaptive radiations, whereas Plethodon is sometimes referred to as a nonadaptive radiation, because many speciation events in this clade appear to have been ecologically equivalent geographical replacements [36]. All three phylogenies clearly show patterns of slowdowns in diversification rate towards the present [22,36,41].
For all three clades, likelihood maximization of the diversitydependent model always leads to nonzero estimates of extinction (table 2). For Dendroica warblers, we find that the DDL+E model provides a substantially better fit to the LTT than a model without extinction (figure 3a and table 2). Under the maximumlikelihood DDL+E model, we predict that this clade reached its carrying capacity 4 Myr ago, after which the speciation and extinction rates have been in equilibrium (figure 3j). We note that the behaviour of our model is similar to Rabosky's [42] model of heritable extinction with pulsed turnover (HEPT), which also yielded a reasonable fit to the Dendroica LTT plot. Both our own and Rabosky's models find that warbler diversity has been constant over the last 4 Myr, with speciation and extinction in equilibrium throughout this period. However, the fiveparameter HEPT model assumes that Dendroica diversity has been constant back to 5 Myr by modelling diversification as a Moran process, whereas our threeparameter model predicts it.
The Plethodon and Heliconius phylogenies are also consistent with diversitydependent speciation, with DDL+E and DDL−E models preferred to the CR model. In both cases, however, the DDL+E model has only marginally higher likelihood than the DDL−E model. In fact, the likelihood profiles for extinction are quite flat across a range of values (figure 3), even though the likelihood surfaces around the optimum are not so flat (errors in table 2): this means that the parameters are highly correlated.
4. Discussion
We have presented a new generalized likelihoodbased inference method that allows estimation of extinction rates from the branching times of extant species when speciation is diversitydependent. This has not previously been possible. The most precise estimation of extinction (and of the other parameters) occurs for LTT plots showing either saturation towards the present or an inverted Sshape (figure 1c and table 1), in which case we expect diversity to be equilibrating in the fossil record. For LTT plots of other shapes, phylogenies do not seem to contain sufficient information to distinguish between various models (figure 1a and table 1), and our predictions for historical diversity are less certain. Nevertheless, even in these cases, our approach brings us much closer to resolving the seeming inconsistency between fossil record and molecular phylogenies: diversitydependence clearly allows for nonzero extinction rates, with the likelihood of the branching times comparable across a wide range of extinction rates. The Akaike weights of our five case studies indicate that the DDL+E outperforms the CR birth–death model in all case studies (for Cetacea when extinction was fixed at the fossil estimate) and the diversitydependent model without extinction in the foraminifera and Dendroica (table 2).
The agreement between our phylogenybased parameter estimates for foraminifera and those obtained by Ezard et al. [15], who used the complete fossil phylogeny, is encouraging. Both studies identify a diversitydependent decline in speciation rate and estimate similar, nonzero, extinction rates. By using the whole tree of extinct as well a extant species, Ezard et al. [15] identified additional complexity in the diversification dynamics, including dependence of dynamics on species' ecology, an influence of the abiotic environment on diversification rates (particularly extinction) and some diversitydependence in the extinction rates. Nonetheless, the fit of the threeparameter DDL+E model to the extant species phylogeny and the prediction of fossil diversity may be difficult to surpass using extant species data alone.
The Dendroica branching times are substantially more likely under the DDL+E model than either the DDL−E or CR models. Note that, while the topology we use for Dendroica is the same as in Rabosky & Lovette [22], our maximumlikelihood parameter estimates differ from theirs because we treat the root as 5 Myr, whereas they scale the tree so that the roottotiplength equals one. The DDL+E model (LL = 16.05) also outperforms the DDX (=DDX − E) model (LL = 13.45) and the SPVAR model (a model where speciation declines through time and extinction is constant, LL = 11.40), which performed best in earlier studies of Rabosky & Lovette [22,23], respectively. Our findings agree with earlier studies that this North American radiation was initially explosive and then slowed down [22,23,35], but our study differs from earlier studies, in that we estimate appreciable extinction rates and a much faster speciation rate at the crown of the tree than that of the DDL−E and SPVAR models (though not the DDX model).
While the foraminifera and Dendroica show an inverted Sshape or a plateau, respectively—in which cases parameter estimates are accurate (table 1)—Cetacea, Heliconius and Plethodon do not show these patterns. We therefore expect the parameter estimates to be uncertain and potentially subject to bias. Indeed, the likelihood surface is relatively flat across a wide range of extinction rates including zero (figure 3e,f) and the errors in the estimates are large (table 2). The Akaike weights of the DDL+E model are not the highest, but still appreciable. From a Bayesian perspective, if there were grounds to favour a low prior probability for zero extinction, the diversitydependent model with extinction would always have outperformed its zeroextinction counterpart, as the posterior probability would recover the prior. However, identifying an appropriate prior distribution for extinction is not straightforward, particularly if extinction is phylogenetically clustered [42].
While our focus has been on modelling extinction (owing to the apparent conflict between phylogenetic estimates and fossil evidence in previous studies), we note that if diversification is indeed diversitydependent then our approach will also provide more accurate estimates of other parameters of diversification. By the same token, if the likelihood surface for extinction is flat, estimates of speciation rates will be subject to similar uncertainty (table 2). We note that our parameter λ_{0} denotes only the initial speciation rate. The speciation rate rapidly decreases as diversity increases, and hence, in most cases (but not in Cetacea), the speciation rate at the present is similar in magnitude to the extinction rate. If we do not allow species other than those at the crown to exist and contribute to diversitydependence at the present, a large initial speciation rate is required to make present speciation and extinction rates similar. As such, our method can also be used to detect the presence of other species at the crown age (see electronic supplementary material), but we have no evidence to suggest that multiple closely related competitors existed at the root of the radiations we consider here.
Various hypotheses other than diversitydependence have been offered to explain macroevolutionary dynamics. In foraminifera, the changing climate appears to play a crucial role; for instance, it may explain the drop in diversity around 34 Myr ago [15]. Likewise, in Cetacea, there is evidence that restructuring of the oceans during the periods 35–31 and 13–4 Myr ago caused temporary increases in speciation rate [33]. We do not deny the influence of external factors on macroevolutionary dynamics, but these external drivers are idiographic [18] and will need to be considered on a casebycase basis. In contrast, we argue that there is now ample and diverse evidence that diversitydependence appears to be a general, nomothetic, phenomenon [18], if not ubiquitous [13–17,19–21] and has a pronounced effect on macroevolutionary dynamics. We suggest that the diversitydependent model with extinction should be preferred to the CR birth–death or pure birth model as a more biologically realistic model for macroevolution. Although the CR birth(–death) model is a more appropriate null model from a statistical perspective because of its greater simplicity, we believe that the generality of diversitydependence requires a more realistic null: models incorporating other (particularly idiographic) mechanisms should be assessed against the (nomothetic) model of diversitydependence. Also the diversitydependent model contains the CR birth and birth–death models as special cases; hence, parameter estimates can inform us whether diversitydependence is strong or weak.
The MEDUSA method [43] has been developed to identify—on the basis of higher level phylogenies and the species richness of clades—the locations of rate transitions in temporally homogeneous constant speciation and extinction rates, with a CR birth–death model as the null. This method could be extended by introducing DDL+E as a more realistic null model and allowing for transitions in each of the three DDL+E parameters. There is no longer a computational obstacle to this use of the DDL+E model.
Diversitydependence need not always be negative. As the clade diversifies, more niches can be constructed, either by external factors (climate), or owing to the mere presence of the species themselves [2,44]. In the former case, this can be incorporated in our framework by enlarging K. In the latter case, the constructed niches do not affect the diversitydependence of the clade under consideration (e.g. arrival of predators [2]) or a key innovation has opened up new opportunities; this can be modelled as a decoupling of diversitydependent dynamics of the innovative subclade from the main clade. This model is explored elsewhere [45].
While our method can estimate nonzero extinction, our estimates for Cetacea, Plethodon and Heliconius are low and likely to underestimate the true rate. This happens because most reconstructed phylogenies do not show the pullofthepresent that is the signature of a high extinction rate [6,7]. When the fact that speciation is not instantaneous is incorporated in the CR birth–death model, the pullofthepresent is diminished, disappears entirely or is even transformed into a concave curve, depending on the time taken for speciation to complete [11]. We conjecture that incorporating protracted speciation in the DDL+E model will result in higher extinction estimates. A likelihood formula is not yet available for protracted speciation in the CR model, let alone in the DDL+E model. Moreover, because our results (figure 1, table 1 and electronic supplementary material, figure S2) indicate that the branching patterns of real phylogenies are equally likely to arise under a wide range of models, the prospects for estimating parameters under other models, which are more complex still, are likely to be limited.
We have studied one of the simplest models of diversitydependence. It may be viewed as a phenomenological description, similar to logistic growth in a population. However, it may also be interpreted more mechanistically, similar to colonization–extinction dynamics in a metapopulation [46]; the underlying assumptions are then that speciation can fill any one of the K′ niches with equal probability and that the ecological dynamics (population growth) of a newly arisen species are much faster than the evolutionary dynamics of diversification. This interpretation opens up a suite of mathematical tools and results of metapopulation ecology, for example, that the predictions for diversity are still good approximations even if species vary in their extinction rates owing to population size fluctuations [47]. Nevertheless, various other diversitydependence models are possible, including ones that do not set a maximum to diversity (e.g. with an exponentially declining speciation rate with diversity). While the latter property seems preferable, it lacks a mechanistic underpinning. We hope that our work will stimulate research in this direction, as we have provided a framework in which computing the likelihood of the phylogenetic branching pattern associated with such models no longer presents a barrier.
Our analyses demonstrate both the potential and limitations of estimating separate speciation and extinction rates from molecular phylogenies. Although the absence of extinct lineages from molecular phylogenies implies that they provide an incomplete picture of macroevolutionary dynamics, we have shown that a more realistic model, which includes diversitydependence and allows for nonzero extinction, returns parameter estimates that are not at odds with the fossil record. Thus, even when fossils are absent, molecular phylogenies can be an informative means of inferring past speciation and extinction dynamics.
Acknowledgements
We thank T. Barraclough, A. Humphreys, T. Quental, D. Rabosky and J. Rosindell for comments on an earlier version of the manuscript, T. Ezard, S. Ho, K. Kozak, J. Mallet, C. Marshall, T. Quental, D. Rabosky and M. Steeman for providing phylogenetic trees or fossil data, and M. Foote for advice. We also thank two anonymous reviewers for constructive comments. The Heliconius and foram photos were provided by J. Mallet and H. Coxall, respectively; the other photos are under Public Domain license on Wikimedia Commons. We thank P. Etienne for photo editing. Financial support for R.S.E. was provided by a VIDI grant from the Netherlands Scientific Organization (NWOALW) and a FRoSpects Short Visit Grant from the European Science Foundation (ESF), for R.S.E. and B.H. by an exchange grant from the Van Gogh Programme, and for T.A. by NERC (grant NE/E015956/1 to A.P. and P.N.P.).
 Received July 12, 2011.
 Accepted September 22, 2011.
 This journal is © 2011 The Royal Society