# How Well Can We Estimate Diversity Dynamics for Clades in Diversity Decline?

How Well Can We Estimate Diversity Dynamics for Clades in Diversity Decline? Abstract The fossil record shows that the vast majority of all species that ever existed are extinct and that most lineages go through an expansion and decline in diversity. However, macroevolutionary analyses based upon molecular phylogenies have difficulty inferring extinction dynamics, raising questions about whether the neontological record can contribute to an understanding of the decline phenomenon. Two recently developed diversification methods for molecular phylogenies (RPANDA and BAMM) incorporate models that theoretically have the capacity to capture decline dynamics by allowing extinction to be higher than speciation. However, the performance of these frameworks over a wide range of decline scenarios has not been studied. Here, we investigate the behavior of these methods under decline scenarios caused by decreasing speciation and increasing extinction through time on simulated trees at fixed intervals over diversity trajectories with expansion and decline phases. We also compared method performance over a comprehensive data set of 214 empirical trees. Our results show that both methods perform equally well when varying speciation rates control decline. When decline was only caused by an increase in extinction rates both methods wrongly assign the variation in net diversification to a drop in speciation, even though the positive gamma values of those trees would suggest otherwise. We also found a tendency for RPANDA to favor increasing extinction and BAMM to favor decreasing speciation as the most common cause of decline in empirical trees. Overall our results shed light on the limitations of both methods, encouraging researchers to carefully interpret the results from diversification studies. BAMM, diversification, diversity decline, extinction, macroevolution, RPANDA, simulations, speciation, vertebrates Accurate assessment of biodiversity changes over geological time is critical to an understanding of the mechanisms that control species richness. One of the most striking patterns of historical biodiversity is the waxing and waning of clades; over 99% of all species that have ever existed are extinct, and the diversity decline phase of clades is a well-recognized phenomenon across the tree of life (Raup 1986; Foote 2007; Liow and Stenseth 2007; Quental and Marshall 2013; Silvestro et al. 2015). Traditionally, paleontologists have assessed extinction and origination using the fossil record to understand the deep time processes that shaped biodiversity. These studies have helped to clarify several extinction-related processes including biodiversity rebounds following mass extinction events (Roopnarine 2006; Stanley 2007), the possibility that diversity might regulate itself (Sepkoski 1978, 1979, 1984; Alroy 1996), that the diversity dynamics of different clades might be interconnected (Liow et al. 2015; Silvestro et al. 2015), and the fact that the diversity decline might be controlled by both a drop in speciation as a rise in extinction (Gilinsky and Bambach 1987; Quental and Marshall 2013). The fossil record provides the most direct way of estimating extinction and speciation through time (Nichols and Pollock 1983; Quental and Marshall 2009, 2010); however, the degree of incompleteness of the record may potentially limit inferences of such rates (Peters 2005; Alroy 2010). Over the last several decades methods that explicitly account for sampling and preservation biases have been developed (e.g., Alroy 1996; Foote 2007; Silvestro et al. 2014; Starrfelt and Liow 2016); however, these methods still require a minimum amount of stratigraphic information. Despite a growing interest among neontologists in incorporating fossil data into phylogenetic studies (Quental and Marshall 2010; Morlon et al. 2011; Slater et al. 2012; Heath et al. 2014), the number of empirical systems that permit the integration or the comparison of molecular phylogenetic data with rich stratigraphic information is small. As a result, there are relatively few examples where the fossil record allows for unambiguous identification of extant clades in decline (e.g., Silvestro et al. 2014). The advent of statistical methods to estimate speciation and extinction rates solely from molecular phylogenies (Nee et al. 1994; Stadler 2013; Morlon 2014; Rabosky 2014) and the unprecedented recent accumulation of molecular phylogenies offers the intriguing possibility that extant phylogenies alone could be used to identify clade decline dynamics. These approaches, ranging from simple constant-rate birth–death models (Nee et al. 1994) to complex Bayesian (Rabosky 2014), and trait-dependent models (FitzJohn 2012), have been heavily used in the past two decades. The use of molecular phylogenies has not only allowed researchers to evaluate the diversification dynamics of several types of organisms with a poor fossil record (e.g., birds Huang and Rabosky 2014; Burin et al. 2016), but also opened up the possibility for broader generalizations due to the paramount quantity of molecular data. However, several authors have cautioned against the ability of birth–death models to properly infer extinction rates (Rabosky 2010, but see Beaulieu and O’Meara 2015 and Rabosky et al. 2015) and hence properly detect diversity decline (Quental and Marshall 2010). With few exceptions, empirical studies using molecular phylogenies suggest either very small extinction rates (e.g., Rabosky 2014; Alencar et al. 2016) or that diversity is either increasing or at best plateauing (e.g., Morlon et al. 2010; Rabosky 2014). One limitation of constant-rate models is that they will always reconstruct positive net diversification rates on molecular phylogenies, since diversity is always increasing over clade history on trees that only include extant taxa. Thus, clade decline patterns should be impossible to detect with these models and/or data sets (Barraclough and Nee 2001; Ricklefs 2007; Quental and Marshall 2010; Quental and Marshall 2011). Two recently developed diversification methods, RPANDA and BAMM, circumvent this limitation by allowing extinction to be higher than speciation (Morlon et al. 2011; Rabosky 2014). RPANDA is a likelihood-based model that allows the fitting and selection of a priori models on selected trees (Morlon et al. 2011; Morlon 2014); this method is implemented in the RPANDA package (Morlon et al. 2016). BAMM (Rabosky 2014) relies on a reversible-jump Markov Chain Monte Carlo (MCMC) algorithm to detect shifts in diversification rates and get average rates for the whole tree, as well as for subclades. A challenge to the estimation of extinction from molecular phylogenies lies in detecting dynamically changing speciation and extinction rates as the diversification process unfolds (Liow et al. 2010; Quental and Marshall 2011). Although RPANDA and BAMM both allow for time-varying models of speciation (and also of extinction in RPANDA), the performance of these methods under different scenarios of changing net diversification rate is relatively poorly understood. RPANDA has been shown to give reasonable performance over a limited range of decline scenarios (Morlon et al. 2011), albeit in a much more restricted parameter space, requiring a more thorough examination of the method in a wider range of scenarios for a better assessment of the true power of the method. It is also worth mentioning that BAMM does not explicitly models varying extinction rate, but given its plausibility, and the fact that most researchers will use BAMM without acknowledging such variation, it is nonetheless important to test its performance on such scenario. It is interesting to note that both RPANDA and BAMM showed contrasting results when analyzing the same empirical tree (Cetacea), and that both empirical analyses show important differences from what is inferred from the fossil record (Quental and Marshall 2010). In their initial study, Morlon et al. (2011) showed that RPANDA recovered a diversity decline for Cetaceans, that is, for the most part driven by a rise in extinction rates, while BAMM suggested expanding diversity (Rabosky 2014). One should note that such comparison is not straightforward given that Morlon et al. (2011), differently than Rabosky (2014), recovered negative diversification rates in Cetacea after splitting up the phylogeny in different subgroups, including paraphyletic ones. The cetacean fossil record on the other hand suggests that a drop in speciation is of primary importance in the diversity decline dynamics not only for cetaceans (Quental and Marshall 2010), but also for many other organisms (Gilinsky and Bambach 1987; Quental and Marshall 2013; Silvestro et al. 2014, 2015). Given that molecular phylogenies allow researchers to only indirectly estimate the contribution of speciation and extinction regimes, it is possible that such methods might have a difficult time in properly estimating their individual changes. A survey on the literature using Google Scholar for papers citing either Rabosky (2014) or Morlon et al. (2011) revealed a surprisingly low number of papers that acknowledge and/or discuss decline in diversity (see Discussion section for a full description of the survey). This low number of studies shows that diversity decline is, in general, overlooked by researchers, and one possible explanation is the inability of the methods (but see Discussion section for a more detailed examination). In contrast, clades declining in diversity comprise a pattern commonly seen in the fossil record. Although such rarity is not direct evidence, given what we know from the fossil record it is unlikely that almost all clades ever studied using molecular phylogenies belong to a dynamics of either expansion or equilibrium diversity (Quental and Marshall 2010). We indeed agree that estimating extinction rates from molecular phylogenies is problematic (Rabosky 2010), and thus the investigation of a scenario where extinction clearly plays an important role, such as diversity decline, deserves full attention, and a thorough investigation. If the current, state-of-the-art birth–death methods can reliably estimate diversification from molecular phylogenies, in particular the ability to distinguish expansion, equilibrium, and decline diversity scenarios, they could allow researchers to scan the tree of life with no taxonomic limitation, opening up the possibility for broad generalizations. More importantly, it would allow the reconciliation and synthesis between the paleontological and neontological records. Therefore, here we present a broad analysis to compare the ability of the two most recent and promising methods on estimating diversification dynamics, with special attention to their ability to detect diversity decline and properly infer the role of speciation and extinction on controlling such dynamics. We address three main questions: 1) What is the performance of RPANDA and BAMM when estimating rates of diversification under different scenarios of decline in diversity? 2) Can these methods discriminate between decline scenarios caused by decreasing speciation rates and increasing extinction rates? 3) Does the ability of each method to properly infer changes in speciation and extinction rates, and their net result, depend on the intensity of the decline? To address these questions, we investigate a range of evolutionary scenarios over a comprehensive parameter space at different points in time in a clade’s history. We also analyzed a comprehensive data set of empirical trees to address how both methods perform in “real world” scenarios, and to investigate if both methods converge on the same answer when the dynamics is likely to be more complex than in our simulations. Materials and Methods Workflow The workflow for the simulation portion of our study consists of four major steps (Fig. 1). First, we simulated 2000 trees for each of two different diversification scenarios. In the first scenario, the decline in diversity comes from varying (decreasing) speciation and constant extinction rates (hereafter called SPvar), and the second scenario consists of constant speciation and varying (increasing) extinction rates (hereafter called EXvar). Throughout this article, we will consider a tree to be in decline in diversity if the estimated extinction rate is higher than the speciation rate at the present (i.e., at the tips of the tree). Second, we perform a “time travel procedure” wherein each simulated tree is trimmed in a predetermined point in time to emulate what a given molecular phylogeny would look like in different phases of the decline process. Third, we fit both methods to the simulated trees, in order to obtain the diversification rates estimates. Finally, we analyze the performance of the models by evaluating the relationship between estimated and simulated rates. Figure 1. View largeDownload slide Schematic workflow of the study. The top panel represents the main steps of tree simulation and rate estimation (with I–IV being the main sections of the protocol described below, and the pruning step representing the removal of the extinct species). The bottom panels show a cartoon example of the diversity trajectory we simulated as well as an example tree showing the time slices used on the study. (80% rising |$=$| 80% of peak diversity in the diversity expansion phase; 80%/50%/20% left |$=$| 80%/50%/20% of peak diversity in the decline phase). Figure 1. View largeDownload slide Schematic workflow of the study. The top panel represents the main steps of tree simulation and rate estimation (with I–IV being the main sections of the protocol described below, and the pruning step representing the removal of the extinct species). The bottom panels show a cartoon example of the diversity trajectory we simulated as well as an example tree showing the time slices used on the study. (80% rising |$=$| 80% of peak diversity in the diversity expansion phase; 80%/50%/20% left |$=$| 80%/50%/20% of peak diversity in the decline phase). I: Simulations — The first step consisted of simulating trees for two different diversification scenarios: the first scenario was characterized by an exponential decline on speciation rates and constant extinction rates through time (SPvar, Table 1), whereas the second scenario consisted of constant speciation rates and increasing extinction rates over time (EXvar, Table 1). We chose the saturating increase instead of an exponential increase for the EXvar scenario so the variation in net diversification between both scenarios would follow the same pattern. In the SPvar scenario, we used two parameters for speciation rates (initial speciation |$\lambda_{0}$| and decaying rate |$\alpha$|⁠) and one parameter for extinction rates |$\mu$|⁠. In the second scenario (EXvar), we used one parameter for speciation rates |$(\lambda)$| and two for extinction rates (initial extinction |$\mu_{0}$| and increase rate |$\beta)$|⁠. Table 1. Models of rate variation from the root to the tips and ranges of the uniform distributions from which each parameter was sampled in the two simulated scenarios (SPvar and Exvar) Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Table 1. Models of rate variation from the root to the tips and ranges of the uniform distributions from which each parameter was sampled in the two simulated scenarios (SPvar and Exvar) Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Based on the range of values shown in Table 1, we simulated 2000 trees for each scenario. For each simulation, we followed the same three steps: 1) the three parameters necessary for the simulations were randomly sampled from uniform distributions with limits described in Table 1. With the combination of the three parameters, we calculated the time expected for the tree to reach its peak diversity, that is, the time when speciation and extinction rates are equal |$(t_{\mathrm{peak}})$|⁠, as well as the expected maximum number of species. 2) We then estimated the time where the tree is expected to have only 20% of its diversity after passing its peak diversity, and designated this the total simulation time. 3) Finally, we ran each individual simulation using their three sampled parameters and the time of 20% after passing its peak diversity. For the simulations we used a script kindly provided by Dr. Hélène Morlon, that is available on github (http://github.com/gburin/bamm_rpanda). Should the simulated tree reach more than 20,000 species at any point in time, or should the tree go extinct before the defined time, the sampled parameters were recorded, discarded, and the Steps 1–3 were repeated to generate all the 2000 viable trees. The final parameter space was also constrained by defining different limits to the uniform distribution used to sample extinction |$(\mu)$| and speciation |$(\lambda)$| rates where distinct between the two scenarios. This was done to ensure a similar change in net diversification over time (Table 1). We also note that in the varying extinction scenario (EXvar) the final values of extinction can be extremely high and a lot higher than many empirical estimates. Although some of those final extinction rate values might seem unrealistic (in the order of 10 species per lineage per million years), from the standpoint of model performance comparison they should not be a problem. For example, one could rescale the process to unfold at 10 times slower by decreasing the rates by a factor of 10. That would make the rate absolute values similar to empirical ones but the relative relationships between the rates would remain similar to the ones presented here. Lastly it should be noted that those extremely high extinction values are only experienced at the end of the simulation process and that extinction is evaluated at different points in time (see II: Time travel, pruning, and summary statistics section below). As a result, the simulations for both scenarios allowed us to explore two different regions of the parameter space. Hence, it is important to note that it is difficult to directly compare the explored parameter spaces between the two scenarios (Fig. 2a, b) because a given rate is either fixed or time-varying depending on the scenario. Figure 2. View largeDownload slide Explored parameter spaces and simulation results. a) Parameter space explored in the SPvar scenario. b) Parameter space explored in the EXvar scenario. c) Distribution of number of tips in the viable trees of both scenarios. d) Distribution of tree depth (in million years) in the viable trees of both scenarios. e) Distribution of the corrected gamma statistic calculated for all the simulated trees from both scenarios. The black points in panels a and b represent the parameter combinations used by Morlon et al. (2011). (Colored figures in online version). Figure 2. View largeDownload slide Explored parameter spaces and simulation results. a) Parameter space explored in the SPvar scenario. b) Parameter space explored in the EXvar scenario. c) Distribution of number of tips in the viable trees of both scenarios. d) Distribution of tree depth (in million years) in the viable trees of both scenarios. e) Distribution of the corrected gamma statistic calculated for all the simulated trees from both scenarios. The black points in panels a and b represent the parameter combinations used by Morlon et al. (2011). (Colored figures in online version). II: Time travel, pruning, and summary statistics — Since the historical information stored in a phylogeny changes through its evolutionary history, our studied assessed model performance at four different phases in the history of the clade. Thus, the following step of the protocol consisted of using the 2000 simulated trees that are in a late stage of decline that corresponds to 20% of its diversity after passing its peak diversity (hereafter called 20left, which refers to the percentage of peak diversity remaining) and time travel back to the expected times when each of these trees would correspond to 50% and 80% of its peak diversity (hereafter called 50left and 80left, respectively, which refers to the percentage of peak diversity left), as well as when the trees would have had 20% less species than at their peak diversity during the expansion phase (hereafter called 80rise, which refers to the percentage of peak diversity attained at a time slice prior to the peak diversity—see Fig. 1 for all time slices). To generate phylogenies that would correspond to those diversities, we estimated the time it would take for trees to have reached 50% and 80% of the peak diversity in the decline phase, and the time for the 80% of peak diversity in the expansion phase. Those time estimates where then used in the function timeSliceTree from the paleotree package in R (Bapst 2012) to slice the 2000 trees in the corresponding periods in the past to generate new trees. The diversity in each time slice was determined by all of the taxa present at that time, even if some of those taxa become extinct in a later time slice. This procedure resulted in the 16,000 phylogenetic trees (2 scenarios * 4 time slices, including the final time * 2000 simulations per scenario) that were used in the final step of the workflow. We characterized all simulated trees at both scenarios and at all time slices using three frequently used features: number of tips, age of the tree, and the gamma statistic (Pybus and Harvey 2000). The gamma statistic is known to be strongly dependent on tree size (number of tips); therefore, we used a version of the gamma value that is standardized in relation to the maximum possible value of gamma (Quental and Marshall 2011) for a tree with the same number of tips, calculated according to equations (1) and (2): \begin{gather} \frac{\gamma}{\gamma_{\mathrm{max}}} \end{gather} (1) \begin{gather} \gamma_{\mathrm{max}}=\frac{\sqrt \frac{1}{12\ast \left(N\,\mathit{tip}-2\right)}}{2} \end{gather} (2) III: Model fitting and parameter estimation — Parameter estimation in RPANDA and BAMM reflect the differences in the underlying statistical frameworks of these methods. RPANDA’s parameter estimation is based on a model selection approach, based on Maximum Likelihood estimation. In this approach, the user provides one or more models (functions) of rate variation through time, and RPANDA estimates the parameter combination that maximizes the likelihood of each model (Morlon et al. 2016). This is valid for both speciation and extinction. In contrast, BAMM is a Bayesian method that uses a reversible-jump MCMC algorithm to explore the posterior distributions of parameters, while it simultaneously tests for the presence of shifts in the diversification regime within the tree and test different diversification models for each regime (Rabosky 2014). For each regime, BAMM estimates the parameters of an exponential increase or decrease in speciation (following the model present in the first cell of Table 1), and the parameters for a constant extinction model. In all cases model and parameter estimation, was restricted to trees with more than two tips (but note that the vast majority of trees had a lot more tips than that; see Fig. 2). Based on the differences between the two methods, we adopted a slightly different approach for each one, which are described in the corresponding sections. RPANDA — To estimate speciation and extinction regimes, we performed model selection using RPANDA by testing all possible pairwise combinations between constant and time-varying speciation and extinction rates, using all the four rate variation models presented in Table 1. We called each of the pairwise combinations as: BOTHcst (both rates constant); SPvar (varying speciation and constant extinction); EXvar (constant speciation and varying extinction); BOTHvar (both rates varying). Thus, for each of the 16,000 trees (which include trees in all time slices) we provided the four model combinations (see Table 1) and selected the best model based on the Akaike’s Information Criterion corrected for finite sample sizes (AICc). We used |$\Delta$|AICc |$>$|2 to pick the best model (Burnham and Anderson 2002). Since one of the tested models was always the exact same model that was used in the simulations, we find this conservative value of 2 to be appropriate. The best model was then used as the best estimate of diversification rate values for each individual tree. RPANDA estimates the values for the parameters at the present and the rate of change for those models where one or both of the rates vary through time. To infer the rates at the origin of the clade, we used the parameter estimates, the corresponded model, and the simulated time for each individual tree of interest. All RPANDA analyses were done in parallel within R environment (R Core Team 2016) using the packages dplyr, foreach, doMC, and RPANDA (Revolution Analytics and Weston 2015a,b; Morlon et al. 2016; Wickham and Francois 2016). BAMM — To estimate the speciation and extinction regimes in BAMM, we calculated the mean speciation and extinction rates at the tips and at the root of the tree, as well as the increase/decay parameter for speciation. The control files containing important information for BAMM such as prior parameters and chain length were set individually for each tree using parameters estimated using the BAMMtools package (Rabosky et al. 2014), and all BAMM analysis were run in parallel using the GNU parallel shell tool (Tange 2011) and the latest BAMM version (2.5; Rabosky 2014). After running BAMM, the resulting files were imported in R and processed to retrieve the rates using the package BAMMtools. We note that the scenario with only varying speciation (SPvar) yields the most direct comparisons between the methods, given that the current version of BAMM does not allow a test for varying extinction rates. Even though the varying extinction scenario (EXvar) is only implemented in RPANDA, we decided to also use BAMM to estimate rates under this scenario to compare the methods, but more importantly to understand the behavior of both methods in a potentially more complicated scenario. This is justified given that it is an evolutionary plausible scenario, and from a given molecular phylogeny there is no a priori justification to exclude such scenario as a potential one. Empirical trees — We fitted both RPANDA and BAMM to 214 empirical trees. The trees used here are the same as the trees used by Lewitus and Morlon (2016), and include trees of all five major vertebrate groups (fishes, amphibians, reptiles, birds, and mammals). For RPANDA, we fit all four models used in the simulation analysis plus a fifth model that represented an exponential increase in both speciation and extinction (hereafter called BOTHexp). This model was chosen because it includes an additional scenario of expanding diversity (which might be plausible in empirical phylogenies). Depending on the parameter estimates it can represent either quasiconstant rates (when the increase parameter is close to 0) or increasing rates (when the increase parameter is different from 0), therefore providing a versatile scenario for testing. For BAMM, we estimated the set of parameters for each tree individually using the BAMMtools package for R. Since there was only one large tree (Muridae, 680 species) for which the “expected number of shifts” parameter should be set different than 1, we used 1 for all trees. All the estimated rates were then summarized following the same protocol as for the simulated trees, the only addition being that for the empirical trees, when using RPANDA, we estimated the values for each parameter via model averaging, in addition to analyzing the results for the best selected model. These estimates consisted in the parameter values in each of the four models averaged by their respective Akaike weights. We do not know which model is the true model that shaped the diversity of each empirical clade, whereas for the simulated data set we knew the true model in each scenario. Therefore, we used the model averaging approach only for the empirical estimates. Lastly, we estimated the correlation between the estimated rates of both methods (speciation, extinction, and net diversification) both at the tips and at the root using linear models to help us better understand the differences between the parameter estimates. IV: Performance analysis — All the parameter estimates for all 16,000 trees were tabulated and used to evaluate the performance of the two methods using Pearson’s correlation tests (Table 5). All parameter estimates shown in main results correspond to the best model for RPANDA, regardless of whether or not it is the single best model. Also, to explore whether recent studies using these methods have uncovered clades in decline, we surveyed the recent literature by analyzing the results from papers that cited either Morlon et al. (2011) and/or Morlon et al. (2016), and Rabosky (2014). Results Simulations Our results indicate that the collection of viable trees allowed us to greatly expand the parameter space explored by Morlon et al. (2011) and explore the total parameter space delimited by the values present in Table 1 in a comprehensive manner (Fig. 2a, b). Our simulations also resulted in trees with similar size and depth distributions for both scenarios (Fig. 2c, d), with trees ranging from 2 to more than 6000 species, and from 3.3 to 100 million years of age. Hence, the simulations of both decline scenarios resulted in a virtually indistinguishable distribution of phylogenies with respect to total time and richness. The branching patterns on the other hand are very different. The standardized gamma statistic distributions of both scenarios have very little overlap (Fig. 2e), being predominantly negative for SPvar and positive for more than 98% of the trees at EXvar. Fitting We present the results of three out of the four time slices (80rise, 80left, and 20left, respectively). The results for the 50left time slice is only presented in the Supplementary Material available on Dryad at http://dx.doi.org/10.5061/dryad.tt10r84 (Supplementary Figs. S1–S4 available on Dryad) because they are qualitatively similar, and intermediate to the 80left (discrete decline) and 20left (pronounced decline). The predominantly positive estimated values of net diversification rates for the 80rise time slice indicate that there is no bias towards falsely detecting decline for both methods in both scenarios (Panels a and d of Supplementary Figs. S1 and S2 available on Dryad), nor any bias in estimating individual rates of speciation and extinction (Figs. 3 and 4a, d, g, and j) during the diversity expansion phase. Hence we focus on the two other time slices, which represent our main interest in diversity decline. Figure 3. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the SPvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). Figure 3. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the SPvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). Figure 4. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the EXvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). Figure 4. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the EXvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). SPvar — RPANDA selects the “true” model for the vast majority of the trees (91% for the 80left time slice and 87.5% for the 20left time slice—Table 2). The good performance of RPANDA on estimating net diversification is due to the nonbiased and accurate estimation of both speciation and extinction rates at present (Fig. 3b, c, g, h, and i and Supplementary Figs. S1, S3, S4, and S5 and Table S5 available on Dryad), with a very small bias for both rates in the smaller trees (Supplementary Figs. S3 and S4 available on Dryad). Extinction rate is frequently estimated to be very close to zero when the BOTHcst model is selected. Thus, RPANDA fails to detect decline for a few trees in both time slices, especially when the best model selected is BOTHcst. RPANDA estimates of net diversification show little bias for both time slices (80left and 20left—Supplementary Figs. S1 and S5, and Table S5 available on Dryad). When the most complex model (BOTHvar) was selected, the estimates of net diversification were highly negative, overestimating the pace of diversity decline for those trees (Supplementary Fig. 1b, c available on Dryad), especially for moderately sized trees. Table 2. Number of trees for which RPANDA choses each of the four tested models as the best model in both scenarios. Top level headings indicate the true simulated model in each scenarios. Numbers in bold show the instance where the correct model was chosen SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) Table 2. Number of trees for which RPANDA choses each of the four tested models as the best model in both scenarios. Top level headings indicate the true simulated model in each scenarios. Numbers in bold show the instance where the correct model was chosen SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) In the SPvar scenario, RPANDA detected a single best model for the majority of trees (Table 3). Additionally, it is worth noticing that even when more than one model were equally likely, the true model was among the best models for the vast majority of trees (Supplementary Table S1 available on Dryad). Similarly, BAMM is able to detect diversity decline for most trees in both time slices (Supplementary Fig. S1e, f available on Dryad), as well as a drop in speciation for at least 96.65% of trees (see Table 4). However, BAMM produced slightly negatively biased estimates for net diversification (Supplementary Fig. S1e, f and S8 available on Dryad). This bias is more evident for small trees (up to 100 tips), and also at the 80left time slice (Supplementary Fig. 8 available on Dryad). Interestingly, both speciation and extinction rates seem to be overestimated for trees up to 1000 tips, and this bias is seen more clearly at the 20left time slice although the number of trees is small and hence biases harder to infer (Fig. 3e, f, k, and l, and Supplementary Figs. S6 and S7 available on Dryad). Surprisingly, our results indicate that this bias seems to be stronger for speciation rates than for extinction rates (Supplementary Figs. S6 and S7 available on Dryad). Nevertheless, the biases in both rates seem to be coupled, that is, the overestimation in both rates for the same tree is proportional, since the net diversification rates are much less biased than each rate alone (Supplementary Fig. S8 available on Dryad). Table 3. Number of trees for which RPANDA choses each number of equally likely models based on the AICc differences (⁠|$\Delta$|AICc|${}>{}$|2) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) Table 3. Number of trees for which RPANDA choses each number of equally likely models based on the AICc differences (⁠|$\Delta$|AICc|${}>{}$|2) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) Table 4. Number of trees for which BAMM infers decrease in speciation rates in each simulated scenario 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) Table 4. Number of trees for which BAMM infers decrease in speciation rates in each simulated scenario 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) The amount of extinction experienced by the clade had a differential effect on methods. For example, the comparison between the 80left and the 20left time slices suggests that RPANDA estimates were largely robust to increased amounts of extinction. The method showed a similar ability to detect negative diversification in both time slices (although the number of trees for which the method chooses the simplest model—BOTHcst—increases at the 20left time slice—Table 2). Conversely, BAMM estimates seemed to benefit from an increased duration of the diversity decline phase and detected diversity decline more frequently at 20left than at 80left time slice (Supplementary Fig. S1e, f available on Dryad). For both methods, the estimated speciation and extinction rates at the root of the trees (reconstructed values for RPANDA and estimated values for BAMM) are similar to the simulated values, regardless of the time slice for both methods (Supplementary Fig. S9 and Table S5 available on Dryad). RPANDA tends to overestimate speciation rates at the root more than BAMM, probably as a consequence of those rate estimates being influenced by the propagation of the error in the estimates of the rate at the present and the decay/increase parameter. This overestimation is greater when the best model selected by RPANDA has variation in extinction rates (BOTHvar). Accordingly, estimated extinction rates at the root are underestimated when the BOTHvar model is selected, being very close to 0. Table 5. Estimated |$R^{2}$| for the linear regression between rate estimates of BAMM and RPANDA for each rate and at each point in time (tip or root). Tip represents estimates at present and “root” at the start of the history of the clades. Values in bold are significant (⁠|$P<0.05$|⁠), and values in both bold and italic are considered poor associations regardless of being significant (⁠|$P<0.05$|⁠). Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 Table 5. Estimated |$R^{2}$| for the linear regression between rate estimates of BAMM and RPANDA for each rate and at each point in time (tip or root). Tip represents estimates at present and “root” at the start of the history of the clades. Values in bold are significant (⁠|$P<0.05$|⁠), and values in both bold and italic are considered poor associations regardless of being significant (⁠|$P<0.05$|⁠). Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 EXvar — The performance of both methods was worse for the varying extinction scenario (EXvar). RPANDA rarely selected the true model as the best model (Table 2); at both the 80left and 20left time slices SPvar was the preferred model for most trees (Table 2), followed by the model with both rates constant. It is also worth noting that the method failed to detect the true simulated model for about 99.6% of the trees at the 80left time slice and for about 98.3% at the 20left time slice. Additionally, in the EXvar scenario RPANDA detects a single best model for much fewer trees than in the SPvar scenario (Table 3). Moreover, for the trees for which one or two models were selected as best (⁠|$\Delta$|AICc |$>$|2), the true (simulated) model was almost never among them (Supplementary Table S1 available on Dryad). Lastly, the estimates under the true model in this scenario (EXvar) were biased towards underestimation for both speciation and extinction at present (Supplementary Fig. S10g–l available on Dryad). This underestimation was more severe for extinction rate, and these biases combined made estimates of net diversification rates at the present to be positive under the true model for the majority of the trees (Supplementary Fig. S10i, l available on Dryad). Overall RPANDA seemed unable to detect varying extinction rates. Although RPANDA detects diversity decline for the EXvar scenario, it did by selecting the wrong model (see Table 2). It is interesting to note that RPANDA showed two distinct biases for net diversification rate estimates: it either estimated net diversification rates to be considerably more negative than the simulated values (regardless of tree size), or it tended to estimate this rate to be a lot less negative, or even positive, for a considerable amount of trees, especially those trees up to 1000 tips (Supplementary Fig. S13 available on Dryad). Nevertheless, the estimates for both speciation and extinction rates seemed nonbiased when we looked at all estimates together (Supplementary Figs. S11 and S12 available on Dryad). This implies that for any given tree RPANDA estimated at least one of the two rates at present to be considerably different from the simulated one. Like RPANDA, BAMM was able to detect negative net diversification rates for most trees, even when BAMM’s implemented model (time-varying speciation with fixed extinction) conflicted with the model used in simulation. However, BAMM estimates of net diversification rates are slightly less biased in all tree sizes (Supplementary Fig. S16 available on Dryad). For BAMM the bias is unidirectional, where estimates of net diversification rates tended to be less negative than the simulated values (Supplementary Fig. S16 available on Dryad). This bias difference is clearer for small trees. Speciation and extinction rates were slightly underestimated at the 20left time slice for trees up to 100 tips (Supplementary Figs. S14 and S15 available on Dryad). Furthermore, the rate of suggesting drops in speciation rates are higher than in the SPvar scenario (at least 74.55% of trees in the two scenarios—Table 4). Overall both methods are more likely to detect diversity decline at the 20left time slice than at the 80left time slice. Interestingly RPANDA detected negative diversification more often than BAMM at both time slices. At 80left time slice, RPANDA detected diversity decline for 1799 trees (89.95%), and BAMM detected decline for 1544 trees (77.2%). At the 20left time slice, RPANDA detected diversity decline for 1840 trees (92%) and BAMM for 1684 trees (84.2%) (Supplementary Fig. S2 available on Dryad). In the EXvar scenario, speciation and extinction rates at present are reasonably well estimated in both methods (represented by the small deviations from the perfect fit line in Fig. 4 and Supplementary Table S5 available on Dryad) although the net diversification rate deviates considerably from the simulated values. However, both speciation and extinction rates at the root do not match the simulated values (Supplementary Figs. S17 and S18 and Table S5 available on Dryad). In fact, both methods infer a considerable drop in speciation and a relatively constant extinction at both time slices for most trees, even though the simulated model had constant speciation rate and variable extinction. Even when RPANDA selected variable extinction as the generating model (145 trees for 80left and 181 trees for 20left—adding up trees for which EXvar and BOTHvar models were selected as best models in Table 2), it rarely provided good estimates of the extinction rate at the root of the tree (Supplementary Fig. 18 available on Dryad); on the other hand, BAMM does not allow for extinction to vary, which forces the method to estimate root extinction to have the same values than in the present. Surveying the literature with Google Scholar we found that, until January 2018, 149 empirical studies cited Rabosky (2014) and applied BAMM, whereas 24 empirical studies cited Morlon et al. (2011) and/or Morlon et al. (2016) and applied RPANDA. Several of these studies, however, provide only speciation rates or shift plots without rate values, rendering impossible to know if diversity declines have been detected. From those studies providing enough information for us to know whether or not they detected decline, 9.4% (14 studies) and 20.8% (5 studies) detected decline using BAMM or RPANDA, respectively (for the complete list of publications, see Supplementary Table S1 available on Dryad). The tabulation using RPANDA, although more favorable in terms of finding diversity decline, represent a smaller sample size. Analysis of Empirical Trees RPANDA indicates a decline in diversity (net diversification smaller than |$-$|0.01) for 31% of the analyzed clades, whereas BAMM detects decline for only 6% of clades from the same data set (see also Supplementary Fig. S19 available on Dryad). For only 3.7% of the trees, a negative net diversification rate was inferred for both methods (see also Supplementary Fig. S19 available on Dryad). Our results also indicate that most empirical trees for which decline was inferred also show negative corrected gamma values (Supplementary Fig. S24 available on Dryad). Speciation rates at present were estimated to be similar between both methods (Fig. 5a and Table 5), especially after removing the trees that RPANDA fitted the BOTHexp model, which produced very high, and sometimes nrealistic, extinction rate estimates (Table 4). One should note that the surplus of trees with negative net diversification rates when using RPANDA comes, from the most part, from those trees that fitted models with increase in both rates, the BOTHexp model (43 out of 67 trees; see also Supplementary Table S2 available on Dryad). Extinction rate estimates at present are not significantly correlated between methods, even after excluding the model BOTHexp (see Fig. 5, see “without BOTHexp” results for “Tip” on Table 5). Although RPANDA suggested a broader variation in extinction rates than BAMM (Fig. 5), many of those rates were estimated to be very small at present (e.g., 113 trees with extinction estimates |$\leqslant$|0.01). BAMM, on the other hand, very rarely estimated extinction rates to be so close to 0 (e.g., 2 trees with extinction estimates |$\leqslant$|0.01). The correlation between the estimates of net diversification at present are, as expected, somewhere in between those of the individual rates (Supplementary Fig. S19 available on Dryad; Table 5). Figure 5. View largeDownload slide Estimated rates for BAMM (⁠|$y$|-axis) and RPANDA (⁠|$x$|-axis) for (a, d) speciation rates at the tips and the root, respectively, and for (b, c, e, f) extinction rates at the tips and the root, respectively. The |$x$|-axis of panels b and e were rescaled from panels c and f, respectively, for better visualization, since there are estimated extinction rates from RPANDA that are greater than 10 events/lineage*myr (up to more than 1500 in the detail of c). The red line represents the perfect correlation (estimates from BAMM and RPANDA are equal), and the blue line represents a linear fit between the two variables. |$R^{2}$| values are shown in Table 5. Speciation rates are similar between the two methods, whereas extinction can be quite distinct between the two methods. (Colored figures in the online version). Figure 5. View largeDownload slide Estimated rates for BAMM (⁠|$y$|-axis) and RPANDA (⁠|$x$|-axis) for (a, d) speciation rates at the tips and the root, respectively, and for (b, c, e, f) extinction rates at the tips and the root, respectively. The |$x$|-axis of panels b and e were rescaled from panels c and f, respectively, for better visualization, since there are estimated extinction rates from RPANDA that are greater than 10 events/lineage*myr (up to more than 1500 in the detail of c). The red line represents the perfect correlation (estimates from BAMM and RPANDA are equal), and the blue line represents a linear fit between the two variables. |$R^{2}$| values are shown in Table 5. Speciation rates are similar between the two methods, whereas extinction can be quite distinct between the two methods. (Colored figures in the online version). Inference of speciation rate dynamics showed a remarkable difference between the two methods. RPANDA suggests that the vast majority of empirical trees (166 out of 214—77%) fitted models where speciation was either constant (123 out of 214—57%) or increasing (43 out of 214—20%) (Fig. 6a). BAMM, on the other hand, suggested that the majority of empirical trees (160 out of 214—75%) showed a relative decrease in speciation (Fig. 6a; most points below zero on the |$y$|-axis). Surprisingly, estimates of extinction rate at the root are more correlated than at the present (Fig. 5e and Table 5). This is partially due to the restriction of initial values related to the particularities of two models (SPvar, BOTHcst) that are implemented in BAMM, and to the fact that three (EXvar, BOTHvar, and BOTHexp) out of the five tested models in RPANDA, are strongly constrained to start with very low (or zero) extinction rates (see Fig. 5). Figure 6. View largeDownload slide Relative speciation rate variation, calculated as [lambda at tips |$-$| lambda at root]/max(lambda at tips, lambda at root) of RPANDA as a function of the same index for BAMM. Positive values indicate an increase in speciation towards the present, whereas negative values indicate decreasing speciation rates through time; 0 indicates constant rates. a) Points are colored according to which model detected decline for each empirical tree; b) Points are colored according to which model was selected as the best model by RPANDA. It is possible to notice that the points that deviate the most from the perfect association (solid line) have models with varying extinction rates as the best model, and that the points that represent a model that is the same as the model used by BAMM fall very close to the red line (green dots on the bottom left part of panel b). (Colored figures in the online version). Figure 6. View largeDownload slide Relative speciation rate variation, calculated as [lambda at tips |$-$| lambda at root]/max(lambda at tips, lambda at root) of RPANDA as a function of the same index for BAMM. Positive values indicate an increase in speciation towards the present, whereas negative values indicate decreasing speciation rates through time; 0 indicates constant rates. a) Points are colored according to which model detected decline for each empirical tree; b) Points are colored according to which model was selected as the best model by RPANDA. It is possible to notice that the points that deviate the most from the perfect association (solid line) have models with varying extinction rates as the best model, and that the points that represent a model that is the same as the model used by BAMM fall very close to the red line (green dots on the bottom left part of panel b). (Colored figures in the online version). Discussion Although the development of new methods (Morlon et al. 2011; Rabosky 2014) might in principle allow estimates of net diversification rates to be negative, the question of whether extinction dynamics can be properly estimated from molecular phylogenies remains contentious (Nee et al. 1994; Purvis 2008; Quental and Marshall 2010; Rabosky 2010; Morlon et al. 2011; Beaulieu and O’Meara 2015; Rabosky 2016; Moore et al. 2016; Rabosky et al. 2017). Our results shed new light on the performance of these methods by considering a wide range of decline scenarios. When looking at simple simulated scenarios, we find that both methods are able to detect decline at the present, but both also fail to characterize the temporal dynamics when extinction rate varies. On the other hand, when looking at empirical phylogenies, we find that estimates from both methods are in general not correlated. By the virtue of knowing the true dynamics, simulation studies are extremely valuable to investigate our ability to detect negative diversification rates based solely on molecular phylogenies, but to date this has been limitedly done to only one of the methods discussed here (Morlon et al. 2011). Simulation Perspective on Our Ability to Infer Negative Diversification Rates Although Morlon et al. (2011) suggested that RPANDA was able to detect diversity decline, their study was limited in several important respects: 1) their parameter space was restrictive when compared to broader ranges of biologically plausible values of speciation and extinction rates; 2) it did not explicitly consider method performance at different intensities of decline (here represented by the different time slices); 3) it only examined performance of a single method. Here, we not only explicitly test the ability of BAMM to detect negative diversification rate for the first time, but we expanded the findings by Morlon et al. (2011) by exploring broader parameter space and different points in time of a given diversification scenario. Unexpectedly, we now show that even when the underlying rate dynamics is wrong (wrong model in the RPANDA analysis; or wrong overall dynamics in BAMM), both methods are able to infer negative net diversification at present, at least for the scenarios simulated here. Moreover, our results indicate that the different phase of a tree’s history (our scenarios at 80left, 50left, and 20left time slices) has only a slight impact on the ability of both models in detecting diversity decline at present in the both simulated scenarios (Supplementary Table S5 available on Dryad). This is surprising given that the phylogenetic overall signature of diversification dynamics (the branching patterns) changes as the diversification processes unfolds and the clade ages, a pattern previously revealed by changes in the gamma statistic (Liow et al. 2010; Quental and Marshall 2011). It should be noted that our results also suggested that inferring reliable negative net diversification rates at the recent past is not a proxy for the reliability of the whole diversification dynamics reconstruction. Can We Properly Infer the Diversification Dynamics of Clades in Decline? Our results suggest that the accuracy of the methods in depicting the temporal dynamics, either by choosing the best model (in the case of RPANDA) or by estimating the rates at both the present and at the root of the trees (in the case of BAMM) varied drastically according to the scenario explored. We clearly show that when extinction was allowed to vary (EXvar scenario), both methods performed worse than when extinction was held constant (SPvar scenario) (Supplementary Table S5 available on Dryad). Given that RPANDA is not limited by preimplemented models (e.g., in our analysis we included varying extinction rate models), its inability to infer varying extinction rates suggests that this poorer performance of both methods might be more related to the limited information inherent to molecular data (e.g., no direct information on extinction) than to the method or implementation used. Under the EXvar scenario, speciation rate estimates at the root (Supplementary Fig. S17 available on Dryad) are significantly higher than the true values for both BAMM and RPANDA. Similarly, both methods suggested a significant drop in speciation rate (Fig. 4 and Supplementary Fig. S17 available on Dryad) when in reality speciation was held constant. Moreover, the prevalence of models chosen by RPANDA with varying speciation in the simulated EXvar scenario indicates that RPANDA detected variation in net diversification but explained it as dynamically changing speciation rate. Given that those methods rely on indirect inference of the past, the amount of information available for both methods becomes thinner as we move towards the root of the tree. In fact, both RPANDA and BAMM showed greater associated error for parameters at the root regardless of the simulated scenario (Supplementary Figs. S9 and S18 available on Dryad), a pattern similar to what was found by Moore et al. (2016). This might not be surprising but we argue that it represents one important aspect of the two methods and (likely) the nature of signal in phylogenetic data sets. Methods based on molecular phylogenies could be seen as binoculars with limited range: they allow researchers to infer important aspects of the recent history of a clade, but it loses power and focus as one tries to look towards the deep past. It is worth noting, the gamma statistic (Pybus and Harvey 2000) for trees under the EXvar scenario are virtually all positive (Fig. 2), as expected by previous simulation studies (Pybus and Harvey 2000; Rabosky and Lovette 2008; Quental and Marshall 2009). Although gamma itself is not informative about whether a clade is in decline, our simulations suggest that the gamma statistic when combined with RPANDA and BAMM results could be helpful in discerning among decline scenarios. For instance, if BAMM and RPANDA suggest a decrease in speciation rates through time, but the gamma value is positive, this might in fact represent an increase in extinction rates instead of the variation in speciation (see Results section, Fig. 2). Comparing the Performance of Both Methods in Empirical Phylogenies Although both methods performed similarly (either well or badly) in our simulation scenarios, it is interesting to note that both inferred very different diversification dynamics when applied to our collection of 214 empirical trees. It is especially interesting to note that both methods only agree with respect to negative diversification rates on about 3.7% of the trees. Of course we do not know the true underlying dynamics for the empirical trees, but the comparison is informative nonetheless because a lack of convergence signals that one (or both) method(s) is providing incorrect estimates. Our empirical analysis also reinforces the idea that inferences at deep time might be more challenging than inferences at the recent present. Both methods converged in similar speciation rate estimates at present but inferred very different speciation rates at the root (Fig. 5). BAMM tends to suggest an overwhelming drop in speciation towards the present, irrespective if the lineages were in diversity decline or not while RPANDA show many different patterns of rate variation. This discrepancy might be driven by the important difference on how flexible the two methods are regarding modeling extinction rates because most discrepant estimates in speciation dynamics come from those cases where RPANDA inferred changes in extinction (Fig. 5b, c). Extinction rate estimates at the presents showed an even stronger discrepancy across the methods (Fig. 5b, c). This is not surprising given that extinction rates can never vary with time in BAMM, but this discrepancy has important consequences when it comes to detecting diversity decline. Depending on the method used, one would have a very different view of how common is diversity decline in the real world (31% vs 6% of clades in diversity decline, for RPANDA and BAMM, respectively). This casts serious doubts on our current ability to estimate extinction rates (Rabosky 2010) and, more importantly, detect diversity decline based on current available methods. Our analysis of empirical data sets indicates a conflict between the speciation dynamics for the set of analyzed trees and results from some prior studies. In our study, results from RPANDA suggest that increasing speciation rates are rather common. On the other hand, both the gamma statistic (McPeek 2008; Phillimore and Price 2008) and a more complex model fitting study (Morlon et al. 2010) suggest that slowdowns are pervasive. It should be noted that the method used in Morlon et al. (2010) did not allow (at least in that particular paper) extinction rates to be higher than speciation rates. RPANDA (Morlon et al. 2011) allows net diversification to be negative, hence our results for the empirical data set lead us to infer that either a drop in speciation might be a common phenomenon and RPANDA misfit speciation dynamics for clades in diversity decline; or that a drop in speciation is not so prevalent as once thought. Lastly, the association between underlying drivers of diversity decline estimated by BAMM and RPANDA and the gamma values might be worth exploring. For BAMM the vast majority of clades in decline show negative gamma values (which would suggest a decrease in speciation rates), and those show a drop in speciation. This association between negative gamma values and a drop in speciation is not as strong for the RPANDA estimates. Whereas most clades in decline show negative gamma values, RPANDA did not assigned scenarios of decrease in speciation for many of them (see empirical summary in Supplementary Material available on Dryad). Those that show inconsistency between gamma values and the role of speciation should be seen with caution (Quental and Marshall 2009, 2011). The literature survey revealed a lack of empirical examples of negative diversification rates in the literature that is quite surprising. Even though this paucity could be related to a lack of studies using such methods, the total number of papers analyzed in our literature survey suggests otherwise. Another potential explanation for this empirical lack of negative diversification rates could be that empirical molecular phylogenies are either very recent radiations, or much more heterogeneous than the ones we typically simulate. If this were the case, the prevalence of molecular phylogenies not showing diversity decline would not be surprising, but rather expected. Empirical trees could for example include a radiating clade that masks a potential diversity decline of a subclade (see Moyne and Neige 2007; or the contrast between Morlon et al. 2011 and Rabosky 2014 for potential examples). We cannot evaluate this for our literature survey, but for the 214 empirical trees analyzed here this does not seem to be the case. Our BAMM analysis did not find any significant shifts for the vast majority of those empirical trees (see Results section). Additionally the empirical trees analyzed here have crown group ages that range from 6.7 myr to 131 myr, with a median value of about 30 myr. It is possible that most of the recent radiations are indeed not in decline, but judging from what we know from the fossil record (e.g., Silvestro et al. 2015; Pires et al. 2017), we would expect many of those older radiations to have had enough time to pass their expanding phase. Although we cannot use the empirical results to evaluate methods performance, their lack of convergence casts serious doubts on our ability to properly estimate the underlying diversification dynamics, especially at deep time. Final Remarks: How Reliable are Estimates of Diversification Dynamics The last decades have witnessed the constant development of new comparative methods (Morlon et al. 2011; FitzJohn 2012; Rabosky 2014), and the turmoil of optimistic and pessimistic views on how well can we infer diversification dynamics solely from molecular phylogenies. This field has recently arrived to its final frontier, the idea that one could detect diversity decline, a scenario where the recent is experiencing more extinction than speciation even though no direct information on extinction is stored in phylogenies. Only one study (Morlon et al. 2011) has explicitly tested the ability of phylogeny-based models on detecting negative diversification rates, and only a few empirical studies have in fact detect this scenario. For the first time to our knowledge, we explicitly tested the ability of BAMM (Rabosky 2014) to detect negative diversification rates. Although our simulation results suggest that these methods might be able to detect negative diversification rates, we also showed that the inferred dynamics (either using BAMM or RPANDA) might be utterly wrong. Additionally, our simulation results and the lack of convergence in the empirical analysis suggested that the empirical pattern of drop in speciation previously thought to be prevalent (Morlon et al. 2010) might in fact result, at least in part, from model “misbehavior.” Hence, we suggest that current methods might have a difficult time on simultaneously inferring deep time dynamics and recent dynamics, especially for clades experiencing diversity decline. Furthermore, based on our results we reinforce that researchers should try to use external sources of information (e.g., fossil record when available) when trying to reconstruct diversification dynamics in deep time. This might be especially relevant if there is external evidence for considerable variation in extinction rates. Although we remain pessimistic about our ability to properly infer the diversification dynamics in deep time, especially for clades experiencing diversity decline, we believe that our results might help map the scenarios where those models are more likely to succeed. Lastly, and perhaps more optimistically, we suggest that current phylogenetic methods might at least properly infer very recent dynamics. Supplementary material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.tt10r84. Funding This work was supported by São Paulo Research Foundation (FAPESP) [grants #2014/03621-9 to G.B., #2012/02038-2 to L.R.V.A. and #2012/04072-3 to T.B.Q.]; and NSF [DEB-0918748 to M.E.A.]. Acknowledgments We would like to thank Dr. Hélène Morlon and Dr. Daniel Rabosky for their helpful comments with RPANDA and BAMM. We further thank Dr. Hélène Morlon for all the assistances throughout the whole development of the article. We also thank Mericien Venzon for the aid in running BAMM. We extend our gratitude to Dr. Eric Lewitus and Dr. Hélène Morlon for kindly providing the empirical trees and all students at LabMeMe (Laboratory of Macroevolution and Macroecology) at the University of São Paulo, Dr. Luke Harmon, and two anonymous reviewers for their valuable comments on the study. References Alencar L.R. , Quental T.B. , Grazziotin F.G. , Alfaro M.L. , Martins M. , Venzon M. , Zaher H. 2016 . Diversification in vipers: phylogenetic relationships, time of divergence and shifts in speciation rates . Mol. Phylogenet. Evol. 105 : 50 – 62 . Google Scholar Crossref Search ADS PubMed Alroy J. 1996 . Constant extinction, constrained diversification, and uncoordinated stasis in North American mammals. Palaeogeogr. Palaeoclimatol. Palaeoecol . 127 ( 1–4 ): 285 – 311 Google Scholar Crossref Search ADS Alroy J. 2010 . Fair sampling of taxonomic richness and unbiased estimation of origination and extinction rates. In: Alroy J., Hunt G., editors. Quantitative methods in paleobiology . The Paleontological Society Papers. Vol. 16. The Paleontological Society . Bethesda, MD, USA. 316pp . Bapst D.W. 2012 . paleotree: an R package for paleontological and phylogenetic analyses of evolution . Methods Ecol. Evol. 3 5 : 803 – 807 . Google Scholar Crossref Search ADS Barraclough T.G. , Nee S. 2001 . Phylogenetics and speciation . Trends Ecol. Evol. 16 7 : 391 – 399 . Google Scholar Crossref Search ADS PubMed Beaulieu J.M. , O’Meara B.C. 2015 . Extinction can be estimated from moderately sized molecular phylogenies . Evolution 69 : 1036 – 1043 . Google Scholar Crossref Search ADS PubMed Burin G., Kissling W.D., Guimarães P.R. , Şekercioğlu Ç.H. , Quental T.B. 2016 . Omnivory in birds is a macroevolutionary sink . Nat. Commun. 7 : 11250 . Google Scholar Crossref Search ADS PubMed Burnham K.P. , Anderson D. 2002 . Model selection and multi-model inference. A practical information-theoretic approach. New York, NY, USA : Springer . ISBN 978-0-387-22456-5. FitzJohn R.G. 2012 . Diversitree: comparative phylogenetic analyses of diversification in R . Methods Ecol. Evol. 3 : 1084 – 1092 . Google Scholar Crossref Search ADS Foote M. 2007 . Symmetric waxing and waning of marine invertebrate genera . Paleobiology 33 4 : 517 – 529 . Google Scholar Crossref Search ADS Gilinsky N.L. Bambach R.K. 1987 . Asymmetrical patterns of origination and extinction in higher taxa . Paleobiology 13 : 427 – 445 . Google Scholar Crossref Search ADS Heath T.A. , Huelsenbeck J.P. , Stadler T. , 2014 . The fossilized birth-death process for coherent calibration of divergent-time estimated . Proc. Natl. Acad. Sci. USA 11 29 : E2957 – 66 . Google Scholar Crossref Search ADS Huang H. , Rabosky D.L. 2014 . Sexual selection and diversification: reexamining the correlation between dichromatism and speciation rate in birds . Am. Nat. 184 : E101 – E114 . Google Scholar Crossref Search ADS PubMed Lewitus E. , Morlon H. , 2016 . Natural constraints to species diversification . PLoS Biol. 14 8 : e1002532 . Liow L.H. , Stenseth N.C. , 2007 . The rise and fall of species: implications for macroevolutionary and macroecological studies . Proc. Biol. Sci. 274 1626 : 2745 – 2752 . Google Scholar Crossref Search ADS PubMed Liow L.H. , Quental T.B. , Marshall C.R. 2010 . When can decreasing diversification rates be detected with molecular phylogenies and the fossil record? Syst . Biol. 59 6 : 646 – 659 . Liow L.H. , Reitan T. , Harnik P.G. 2015 . Ecological interactions on macroevolutionary time scales: clams and brachiopods are more than ships that pass in the night. Ecol. Lett. 18 10 : 1030 – 1039 . Google Scholar Crossref Search ADS PubMed McPeek M.A. 2008 . The ecological dynamics of clade diversification and community assembly . Am. Nat. 172 6 : E270 – 84 . Google Scholar Crossref Search ADS PubMed Moore B.R. , Höhna S. , May M.R. , Rannala B. , Huelsenbeck J.P. 2016 . Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures . Proc. Natl. Acad. Sci. USA 113 34 : 9569 – 9574 . Google Scholar Crossref Search ADS Morlon H. 2014 . Phylogenetic approaches for studying diversification . Ecol. Lett. 17 : 508 – 525 . Google Scholar Crossref Search ADS PubMed Morlon H. , Parsons T.L. , Plotkin J.B. 2011 . Reconciling molecular phylogenies with the fossil record . Proc. Natl. Acad. Sci. USA 108 : 16327 – 16332 . Google Scholar Crossref Search ADS Morlon H. , Potts M.D. , Plotkin J.B. 2010 . Inferring the dynamics of diversification: a coalescent approach. PLoS Biol. 8 9 : e1000493 . Morlon H. , Lewitus E. , Condamine F.L. , Manceau M. , Clavel J. , Drury J. 2016 . RPANDA: an R package for macroevolutionary analyses on phylogenetic trees . Methods Ecol. Evol. 7 : 589 – 597 . Google Scholar Crossref Search ADS Moyne S. , Neige P. 2007 . The space-time relationship of taxonomic diversity and morphological disparity in the Middle Jurassic ammonite radiation . Palaeogeogr. Palaeoclimatol. Palaeoecol. 248 : 82 – 95 . Google Scholar Crossref Search ADS Nee S. , May R.M. , Harvey P.H. 1994 . The reconstructed evolutionary process . Philos. Trans. R. Soc. Lond. B Biol. Sci. 344 : 305 – 311 . Google Scholar Crossref Search ADS PubMed Nichols J.D. Pollock K.H. 1983 . Estimating taxonomic diversity, extinction rates, and speciation rates from fossil data using capture-recapture models . Paleobiology 9 2 : 150 – 163 . Google Scholar Crossref Search ADS Peters S.E. 2005 . Geologic constraints on the macroevolutionary history of marine animals. Proc. Natl. Acad. Sci. USA 102 35 : 12326 – 12331 . Google Scholar Crossref Search ADS Phillimore A.B. , Price T.D. 2008 . Density-dependent cladogenesis in birds. PLoS Biol. 6 3 : e71 . Google Scholar Crossref Search ADS PubMed Pires M.M. , Silvestro D. , Quental T.B. 2017 . Interactions within and between clades shaped the diversification of terrestrial carnivores . Evolution 71 7 : 1855 – 1864 . Google Scholar Crossref Search ADS PubMed Purvis A. 2008 . Phylogenetic approaches to the study of extinction . Annu. Rev. Ecol. Evol. Syst. 39 : 301 – 319 . Google Scholar Crossref Search ADS Pybus O.G. , Harvey P.H. 2000 . Testing macro-evolutionary models using incomplete molecular phylogenies . Proc. Biol. Sci. 267 : 2267 – 2272 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2009 . Extinction during evolutionary radiations: reconciling the fossil record with molecular phylogenies . Evolution 63 12 : 3158 – 3167 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2010 . Diversity dynamics: molecular phylogenies need the fossil record . Trends Ecol. Evol. (Amst.) 25 8 : 434 – 441 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2011 . The molecular phylogenetic signature of clades in decline . PLoS One 6 : e25780 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2013 . How the red queen drives terrestrial mammals to extinction . Science 341 : 290 – 292 . Google Scholar Crossref Search ADS PubMed R Core Team 2016 . R: a language and environment for statistical computing. Vienna, Austria : R Foundation for Statistical Computing. Available from: https://www.R-project.org. Rabosky D.L. 2010 . Extinction rates should not be estimated from molecular phylogenies . Evolution 64 : 1816 – 1824 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. 2014 . Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS One 9 2 : e89543 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. 2016 . Challenges in the estimation of extinction from molecular phylogenies: a response to Beaulieu and O’Meara . Evolution 70 1 : 218 – 228 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. , Lovette I.J. 2008 . Density-dependent diversification in North American wood warblers . Proc. Biol. Sci. 275 1649 : 2363 – 2371 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. , Grundler M. , Anderson C. , Title P. , Shi J.J. , Brown J.W. , Huang H. , Larson J.G. 2014 . BAMMtools: an R package for the analysis of evolutionary dynamics on phylogenetic trees . Methods Ecol. Evol. 5 7 : 701 – 707 . Google Scholar Crossref Search ADS Rabosky D.L. , Title P.O. , Huang H. 2015 . Minimal effects of latitude on present-day speciation rates in new world birds . Proc. Biol. Sci. 282 : 20142889 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. , Mitchell J.S. , Chang J. 2017 . Is BAMM flawed? Theoretical and practical concerns in the analysis of multi-rate diversification models. Syst Biol. 66 : 477 – 498 . Google Scholar Crossref Search ADS PubMed Raup D.M. 1986 . Biological extinction in earth history . Science 231 : 528 – 1533 . Revolution Analytics , Weston S. , 2015a . foreach: ProvidesForeach Looping Construct for R. R package version 1.4.3. Available from: https://CRAN.R-project.org/package=foreach. Revolution Analytics, Weston S., 2015b . doMC: foreach parallel adaptor for ‘parallel’. R package version 1.3.4. Available from: https://CRAN.R-project.org/package=doMC. Ricklefs R.E. 2007 . Estimating diversification rates from phylogenetic information . Trends Ecol. Evol. 22 11 : 601 – 610 . Google Scholar Crossref Search ADS PubMed Roopnarine P.D. 2006 . Extinction cascades and catastrophe in ancient food webs . Paleobiology 32 1 : 1 – 19 . Google Scholar Crossref Search ADS Sakamoto M. , Benton M.J. , Venditti C. 2016 . Dinosaurs in decline tens of millions of years before their final extinction . Proc. Natl. Acad. Sci. USA 113 18 : 5036 – 5040 . Google Scholar Crossref Search ADS Sepkoski J.J. Jr. 1978 . A kinetic model of phanerozoic taxonomic diversity I . Analysis of marine orders. Paleobiology 4 3 : 223 – 251 . Google Scholar Crossref Search ADS Sepkoski J.J. Jr. 1979 . A kinetic model of phanerozoic taxonomic diversity II . Early phanerozoic families and multiple equilibria. Paleobiology 5 3 : 222 – 251 . Sepkoski J.J. Jr , 1984 . A kinetic model of phanerozoic taxonomic diversity . III. Post-paleozoic families and mass extinctions. Paleobiology 10 2 : 246 – 267 . Silvestro D. , Schnitzler J. , Liow L.H. , Antonelli A. , Salamin N. 2014 . Bayesian estimation of speciation and extinction from incomplete fossil occurrence data . Syst. Biol. 63 : 349 – 367 . Google Scholar Crossref Search ADS PubMed Silvestro D. , Antonelli A. , Salamin N. , Quental T.B. 2015 . The role of clade competition in the diversification of North American canids . Proc. Natl. Acad. Sci. USA 112 : 8684 – 8689 . Google Scholar Crossref Search ADS Slater G.J. , Harmon L.J. , Alfaro M.E. 2012 . Integrating fossils with molecular phylogenies improves inference of trait evolution . Evolution 66 : 3931 – 3944 . Google Scholar Crossref Search ADS PubMed Stadler T. 2013 . Recovering speciation and extinction dynamics based on phylogenies . J. Evol. Biol. 26 6 : 1203 – 1219 . Google Scholar Crossref Search ADS PubMed Stanley S.M. 2007 . An analysis of the history of marine animal diversity Paleobiology 33 ( sp6 ): 1 – 55 . Google Scholar Crossref Search ADS Starrfelt J. , Liow L.H. 2016 . How many dinosaur species were there? Fossil bias and true richness estimated using a Poisson sampling model. Philos. Trans. R. Soc. Lond. B Biol. Sci. 371 1691 : 20150219 . Tange O. 2011 . GNU parallel—the command-line power tool . login: USENIX Mag. 1 36 : 42 – 47 . Wickham H. , Francois R. 2016 . dplyr: a grammar of data manipulation. R package version 0.5.0 . Available from: https://CRAN.R-project.org/package=dplyr. © The Author(s) 2018. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Systematic Biology Oxford University Press

# How Well Can We Estimate Diversity Dynamics for Clades in Diversity Decline?

Systematic Biology, Volume 68 (1) – Jan 1, 2019
16 pages

Publisher
Oxford University Press
ISSN
1063-5157
eISSN
1076-836X
D.O.I.
10.1093/sysbio/syy037
Publisher site
See Article on Publisher Site

### Abstract

Abstract The fossil record shows that the vast majority of all species that ever existed are extinct and that most lineages go through an expansion and decline in diversity. However, macroevolutionary analyses based upon molecular phylogenies have difficulty inferring extinction dynamics, raising questions about whether the neontological record can contribute to an understanding of the decline phenomenon. Two recently developed diversification methods for molecular phylogenies (RPANDA and BAMM) incorporate models that theoretically have the capacity to capture decline dynamics by allowing extinction to be higher than speciation. However, the performance of these frameworks over a wide range of decline scenarios has not been studied. Here, we investigate the behavior of these methods under decline scenarios caused by decreasing speciation and increasing extinction through time on simulated trees at fixed intervals over diversity trajectories with expansion and decline phases. We also compared method performance over a comprehensive data set of 214 empirical trees. Our results show that both methods perform equally well when varying speciation rates control decline. When decline was only caused by an increase in extinction rates both methods wrongly assign the variation in net diversification to a drop in speciation, even though the positive gamma values of those trees would suggest otherwise. We also found a tendency for RPANDA to favor increasing extinction and BAMM to favor decreasing speciation as the most common cause of decline in empirical trees. Overall our results shed light on the limitations of both methods, encouraging researchers to carefully interpret the results from diversification studies. BAMM, diversification, diversity decline, extinction, macroevolution, RPANDA, simulations, speciation, vertebrates Accurate assessment of biodiversity changes over geological time is critical to an understanding of the mechanisms that control species richness. One of the most striking patterns of historical biodiversity is the waxing and waning of clades; over 99% of all species that have ever existed are extinct, and the diversity decline phase of clades is a well-recognized phenomenon across the tree of life (Raup 1986; Foote 2007; Liow and Stenseth 2007; Quental and Marshall 2013; Silvestro et al. 2015). Traditionally, paleontologists have assessed extinction and origination using the fossil record to understand the deep time processes that shaped biodiversity. These studies have helped to clarify several extinction-related processes including biodiversity rebounds following mass extinction events (Roopnarine 2006; Stanley 2007), the possibility that diversity might regulate itself (Sepkoski 1978, 1979, 1984; Alroy 1996), that the diversity dynamics of different clades might be interconnected (Liow et al. 2015; Silvestro et al. 2015), and the fact that the diversity decline might be controlled by both a drop in speciation as a rise in extinction (Gilinsky and Bambach 1987; Quental and Marshall 2013). The fossil record provides the most direct way of estimating extinction and speciation through time (Nichols and Pollock 1983; Quental and Marshall 2009, 2010); however, the degree of incompleteness of the record may potentially limit inferences of such rates (Peters 2005; Alroy 2010). Over the last several decades methods that explicitly account for sampling and preservation biases have been developed (e.g., Alroy 1996; Foote 2007; Silvestro et al. 2014; Starrfelt and Liow 2016); however, these methods still require a minimum amount of stratigraphic information. Despite a growing interest among neontologists in incorporating fossil data into phylogenetic studies (Quental and Marshall 2010; Morlon et al. 2011; Slater et al. 2012; Heath et al. 2014), the number of empirical systems that permit the integration or the comparison of molecular phylogenetic data with rich stratigraphic information is small. As a result, there are relatively few examples where the fossil record allows for unambiguous identification of extant clades in decline (e.g., Silvestro et al. 2014). The advent of statistical methods to estimate speciation and extinction rates solely from molecular phylogenies (Nee et al. 1994; Stadler 2013; Morlon 2014; Rabosky 2014) and the unprecedented recent accumulation of molecular phylogenies offers the intriguing possibility that extant phylogenies alone could be used to identify clade decline dynamics. These approaches, ranging from simple constant-rate birth–death models (Nee et al. 1994) to complex Bayesian (Rabosky 2014), and trait-dependent models (FitzJohn 2012), have been heavily used in the past two decades. The use of molecular phylogenies has not only allowed researchers to evaluate the diversification dynamics of several types of organisms with a poor fossil record (e.g., birds Huang and Rabosky 2014; Burin et al. 2016), but also opened up the possibility for broader generalizations due to the paramount quantity of molecular data. However, several authors have cautioned against the ability of birth–death models to properly infer extinction rates (Rabosky 2010, but see Beaulieu and O’Meara 2015 and Rabosky et al. 2015) and hence properly detect diversity decline (Quental and Marshall 2010). With few exceptions, empirical studies using molecular phylogenies suggest either very small extinction rates (e.g., Rabosky 2014; Alencar et al. 2016) or that diversity is either increasing or at best plateauing (e.g., Morlon et al. 2010; Rabosky 2014). One limitation of constant-rate models is that they will always reconstruct positive net diversification rates on molecular phylogenies, since diversity is always increasing over clade history on trees that only include extant taxa. Thus, clade decline patterns should be impossible to detect with these models and/or data sets (Barraclough and Nee 2001; Ricklefs 2007; Quental and Marshall 2010; Quental and Marshall 2011). Two recently developed diversification methods, RPANDA and BAMM, circumvent this limitation by allowing extinction to be higher than speciation (Morlon et al. 2011; Rabosky 2014). RPANDA is a likelihood-based model that allows the fitting and selection of a priori models on selected trees (Morlon et al. 2011; Morlon 2014); this method is implemented in the RPANDA package (Morlon et al. 2016). BAMM (Rabosky 2014) relies on a reversible-jump Markov Chain Monte Carlo (MCMC) algorithm to detect shifts in diversification rates and get average rates for the whole tree, as well as for subclades. A challenge to the estimation of extinction from molecular phylogenies lies in detecting dynamically changing speciation and extinction rates as the diversification process unfolds (Liow et al. 2010; Quental and Marshall 2011). Although RPANDA and BAMM both allow for time-varying models of speciation (and also of extinction in RPANDA), the performance of these methods under different scenarios of changing net diversification rate is relatively poorly understood. RPANDA has been shown to give reasonable performance over a limited range of decline scenarios (Morlon et al. 2011), albeit in a much more restricted parameter space, requiring a more thorough examination of the method in a wider range of scenarios for a better assessment of the true power of the method. It is also worth mentioning that BAMM does not explicitly models varying extinction rate, but given its plausibility, and the fact that most researchers will use BAMM without acknowledging such variation, it is nonetheless important to test its performance on such scenario. It is interesting to note that both RPANDA and BAMM showed contrasting results when analyzing the same empirical tree (Cetacea), and that both empirical analyses show important differences from what is inferred from the fossil record (Quental and Marshall 2010). In their initial study, Morlon et al. (2011) showed that RPANDA recovered a diversity decline for Cetaceans, that is, for the most part driven by a rise in extinction rates, while BAMM suggested expanding diversity (Rabosky 2014). One should note that such comparison is not straightforward given that Morlon et al. (2011), differently than Rabosky (2014), recovered negative diversification rates in Cetacea after splitting up the phylogeny in different subgroups, including paraphyletic ones. The cetacean fossil record on the other hand suggests that a drop in speciation is of primary importance in the diversity decline dynamics not only for cetaceans (Quental and Marshall 2010), but also for many other organisms (Gilinsky and Bambach 1987; Quental and Marshall 2013; Silvestro et al. 2014, 2015). Given that molecular phylogenies allow researchers to only indirectly estimate the contribution of speciation and extinction regimes, it is possible that such methods might have a difficult time in properly estimating their individual changes. A survey on the literature using Google Scholar for papers citing either Rabosky (2014) or Morlon et al. (2011) revealed a surprisingly low number of papers that acknowledge and/or discuss decline in diversity (see Discussion section for a full description of the survey). This low number of studies shows that diversity decline is, in general, overlooked by researchers, and one possible explanation is the inability of the methods (but see Discussion section for a more detailed examination). In contrast, clades declining in diversity comprise a pattern commonly seen in the fossil record. Although such rarity is not direct evidence, given what we know from the fossil record it is unlikely that almost all clades ever studied using molecular phylogenies belong to a dynamics of either expansion or equilibrium diversity (Quental and Marshall 2010). We indeed agree that estimating extinction rates from molecular phylogenies is problematic (Rabosky 2010), and thus the investigation of a scenario where extinction clearly plays an important role, such as diversity decline, deserves full attention, and a thorough investigation. If the current, state-of-the-art birth–death methods can reliably estimate diversification from molecular phylogenies, in particular the ability to distinguish expansion, equilibrium, and decline diversity scenarios, they could allow researchers to scan the tree of life with no taxonomic limitation, opening up the possibility for broad generalizations. More importantly, it would allow the reconciliation and synthesis between the paleontological and neontological records. Therefore, here we present a broad analysis to compare the ability of the two most recent and promising methods on estimating diversification dynamics, with special attention to their ability to detect diversity decline and properly infer the role of speciation and extinction on controlling such dynamics. We address three main questions: 1) What is the performance of RPANDA and BAMM when estimating rates of diversification under different scenarios of decline in diversity? 2) Can these methods discriminate between decline scenarios caused by decreasing speciation rates and increasing extinction rates? 3) Does the ability of each method to properly infer changes in speciation and extinction rates, and their net result, depend on the intensity of the decline? To address these questions, we investigate a range of evolutionary scenarios over a comprehensive parameter space at different points in time in a clade’s history. We also analyzed a comprehensive data set of empirical trees to address how both methods perform in “real world” scenarios, and to investigate if both methods converge on the same answer when the dynamics is likely to be more complex than in our simulations. Materials and Methods Workflow The workflow for the simulation portion of our study consists of four major steps (Fig. 1). First, we simulated 2000 trees for each of two different diversification scenarios. In the first scenario, the decline in diversity comes from varying (decreasing) speciation and constant extinction rates (hereafter called SPvar), and the second scenario consists of constant speciation and varying (increasing) extinction rates (hereafter called EXvar). Throughout this article, we will consider a tree to be in decline in diversity if the estimated extinction rate is higher than the speciation rate at the present (i.e., at the tips of the tree). Second, we perform a “time travel procedure” wherein each simulated tree is trimmed in a predetermined point in time to emulate what a given molecular phylogeny would look like in different phases of the decline process. Third, we fit both methods to the simulated trees, in order to obtain the diversification rates estimates. Finally, we analyze the performance of the models by evaluating the relationship between estimated and simulated rates. Figure 1. View largeDownload slide Schematic workflow of the study. The top panel represents the main steps of tree simulation and rate estimation (with I–IV being the main sections of the protocol described below, and the pruning step representing the removal of the extinct species). The bottom panels show a cartoon example of the diversity trajectory we simulated as well as an example tree showing the time slices used on the study. (80% rising |$=$| 80% of peak diversity in the diversity expansion phase; 80%/50%/20% left |$=$| 80%/50%/20% of peak diversity in the decline phase). Figure 1. View largeDownload slide Schematic workflow of the study. The top panel represents the main steps of tree simulation and rate estimation (with I–IV being the main sections of the protocol described below, and the pruning step representing the removal of the extinct species). The bottom panels show a cartoon example of the diversity trajectory we simulated as well as an example tree showing the time slices used on the study. (80% rising |$=$| 80% of peak diversity in the diversity expansion phase; 80%/50%/20% left |$=$| 80%/50%/20% of peak diversity in the decline phase). I: Simulations — The first step consisted of simulating trees for two different diversification scenarios: the first scenario was characterized by an exponential decline on speciation rates and constant extinction rates through time (SPvar, Table 1), whereas the second scenario consisted of constant speciation rates and increasing extinction rates over time (EXvar, Table 1). We chose the saturating increase instead of an exponential increase for the EXvar scenario so the variation in net diversification between both scenarios would follow the same pattern. In the SPvar scenario, we used two parameters for speciation rates (initial speciation |$\lambda_{0}$| and decaying rate |$\alpha$|⁠) and one parameter for extinction rates |$\mu$|⁠. In the second scenario (EXvar), we used one parameter for speciation rates |$(\lambda)$| and two for extinction rates (initial extinction |$\mu_{0}$| and increase rate |$\beta)$|⁠. Table 1. Models of rate variation from the root to the tips and ranges of the uniform distributions from which each parameter was sampled in the two simulated scenarios (SPvar and Exvar) Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Table 1. Models of rate variation from the root to the tips and ranges of the uniform distributions from which each parameter was sampled in the two simulated scenarios (SPvar and Exvar) Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Rate variation Parameter range Speciation Extinction |$\lambda$| |$\alpha$| |$\mu$| |$\beta$| SPvar |$\lambda (t)=\lambda_{0}e^{-\alpha t}$| |$\mu (t)=\mu$| [0 - 10] (initial) [0 - 1] [0 - 2] (constant) NA EXvar |$\lambda (t)=\lambda$| |$\mu (t)=\mu_{f}-({\mu_{f}e}^{-\beta t})$| [0 - 10] (constant) NA [0 - 10] (final) [0 - 1] Based on the range of values shown in Table 1, we simulated 2000 trees for each scenario. For each simulation, we followed the same three steps: 1) the three parameters necessary for the simulations were randomly sampled from uniform distributions with limits described in Table 1. With the combination of the three parameters, we calculated the time expected for the tree to reach its peak diversity, that is, the time when speciation and extinction rates are equal |$(t_{\mathrm{peak}})$|⁠, as well as the expected maximum number of species. 2) We then estimated the time where the tree is expected to have only 20% of its diversity after passing its peak diversity, and designated this the total simulation time. 3) Finally, we ran each individual simulation using their three sampled parameters and the time of 20% after passing its peak diversity. For the simulations we used a script kindly provided by Dr. Hélène Morlon, that is available on github (http://github.com/gburin/bamm_rpanda). Should the simulated tree reach more than 20,000 species at any point in time, or should the tree go extinct before the defined time, the sampled parameters were recorded, discarded, and the Steps 1–3 were repeated to generate all the 2000 viable trees. The final parameter space was also constrained by defining different limits to the uniform distribution used to sample extinction |$(\mu)$| and speciation |$(\lambda)$| rates where distinct between the two scenarios. This was done to ensure a similar change in net diversification over time (Table 1). We also note that in the varying extinction scenario (EXvar) the final values of extinction can be extremely high and a lot higher than many empirical estimates. Although some of those final extinction rate values might seem unrealistic (in the order of 10 species per lineage per million years), from the standpoint of model performance comparison they should not be a problem. For example, one could rescale the process to unfold at 10 times slower by decreasing the rates by a factor of 10. That would make the rate absolute values similar to empirical ones but the relative relationships between the rates would remain similar to the ones presented here. Lastly it should be noted that those extremely high extinction values are only experienced at the end of the simulation process and that extinction is evaluated at different points in time (see II: Time travel, pruning, and summary statistics section below). As a result, the simulations for both scenarios allowed us to explore two different regions of the parameter space. Hence, it is important to note that it is difficult to directly compare the explored parameter spaces between the two scenarios (Fig. 2a, b) because a given rate is either fixed or time-varying depending on the scenario. Figure 2. View largeDownload slide Explored parameter spaces and simulation results. a) Parameter space explored in the SPvar scenario. b) Parameter space explored in the EXvar scenario. c) Distribution of number of tips in the viable trees of both scenarios. d) Distribution of tree depth (in million years) in the viable trees of both scenarios. e) Distribution of the corrected gamma statistic calculated for all the simulated trees from both scenarios. The black points in panels a and b represent the parameter combinations used by Morlon et al. (2011). (Colored figures in online version). Figure 2. View largeDownload slide Explored parameter spaces and simulation results. a) Parameter space explored in the SPvar scenario. b) Parameter space explored in the EXvar scenario. c) Distribution of number of tips in the viable trees of both scenarios. d) Distribution of tree depth (in million years) in the viable trees of both scenarios. e) Distribution of the corrected gamma statistic calculated for all the simulated trees from both scenarios. The black points in panels a and b represent the parameter combinations used by Morlon et al. (2011). (Colored figures in online version). II: Time travel, pruning, and summary statistics — Since the historical information stored in a phylogeny changes through its evolutionary history, our studied assessed model performance at four different phases in the history of the clade. Thus, the following step of the protocol consisted of using the 2000 simulated trees that are in a late stage of decline that corresponds to 20% of its diversity after passing its peak diversity (hereafter called 20left, which refers to the percentage of peak diversity remaining) and time travel back to the expected times when each of these trees would correspond to 50% and 80% of its peak diversity (hereafter called 50left and 80left, respectively, which refers to the percentage of peak diversity left), as well as when the trees would have had 20% less species than at their peak diversity during the expansion phase (hereafter called 80rise, which refers to the percentage of peak diversity attained at a time slice prior to the peak diversity—see Fig. 1 for all time slices). To generate phylogenies that would correspond to those diversities, we estimated the time it would take for trees to have reached 50% and 80% of the peak diversity in the decline phase, and the time for the 80% of peak diversity in the expansion phase. Those time estimates where then used in the function timeSliceTree from the paleotree package in R (Bapst 2012) to slice the 2000 trees in the corresponding periods in the past to generate new trees. The diversity in each time slice was determined by all of the taxa present at that time, even if some of those taxa become extinct in a later time slice. This procedure resulted in the 16,000 phylogenetic trees (2 scenarios * 4 time slices, including the final time * 2000 simulations per scenario) that were used in the final step of the workflow. We characterized all simulated trees at both scenarios and at all time slices using three frequently used features: number of tips, age of the tree, and the gamma statistic (Pybus and Harvey 2000). The gamma statistic is known to be strongly dependent on tree size (number of tips); therefore, we used a version of the gamma value that is standardized in relation to the maximum possible value of gamma (Quental and Marshall 2011) for a tree with the same number of tips, calculated according to equations (1) and (2): \begin{gather} \frac{\gamma}{\gamma_{\mathrm{max}}} \end{gather} (1) \begin{gather} \gamma_{\mathrm{max}}=\frac{\sqrt \frac{1}{12\ast \left(N\,\mathit{tip}-2\right)}}{2} \end{gather} (2) III: Model fitting and parameter estimation — Parameter estimation in RPANDA and BAMM reflect the differences in the underlying statistical frameworks of these methods. RPANDA’s parameter estimation is based on a model selection approach, based on Maximum Likelihood estimation. In this approach, the user provides one or more models (functions) of rate variation through time, and RPANDA estimates the parameter combination that maximizes the likelihood of each model (Morlon et al. 2016). This is valid for both speciation and extinction. In contrast, BAMM is a Bayesian method that uses a reversible-jump MCMC algorithm to explore the posterior distributions of parameters, while it simultaneously tests for the presence of shifts in the diversification regime within the tree and test different diversification models for each regime (Rabosky 2014). For each regime, BAMM estimates the parameters of an exponential increase or decrease in speciation (following the model present in the first cell of Table 1), and the parameters for a constant extinction model. In all cases model and parameter estimation, was restricted to trees with more than two tips (but note that the vast majority of trees had a lot more tips than that; see Fig. 2). Based on the differences between the two methods, we adopted a slightly different approach for each one, which are described in the corresponding sections. RPANDA — To estimate speciation and extinction regimes, we performed model selection using RPANDA by testing all possible pairwise combinations between constant and time-varying speciation and extinction rates, using all the four rate variation models presented in Table 1. We called each of the pairwise combinations as: BOTHcst (both rates constant); SPvar (varying speciation and constant extinction); EXvar (constant speciation and varying extinction); BOTHvar (both rates varying). Thus, for each of the 16,000 trees (which include trees in all time slices) we provided the four model combinations (see Table 1) and selected the best model based on the Akaike’s Information Criterion corrected for finite sample sizes (AICc). We used |$\Delta$|AICc |$>$|2 to pick the best model (Burnham and Anderson 2002). Since one of the tested models was always the exact same model that was used in the simulations, we find this conservative value of 2 to be appropriate. The best model was then used as the best estimate of diversification rate values for each individual tree. RPANDA estimates the values for the parameters at the present and the rate of change for those models where one or both of the rates vary through time. To infer the rates at the origin of the clade, we used the parameter estimates, the corresponded model, and the simulated time for each individual tree of interest. All RPANDA analyses were done in parallel within R environment (R Core Team 2016) using the packages dplyr, foreach, doMC, and RPANDA (Revolution Analytics and Weston 2015a,b; Morlon et al. 2016; Wickham and Francois 2016). BAMM — To estimate the speciation and extinction regimes in BAMM, we calculated the mean speciation and extinction rates at the tips and at the root of the tree, as well as the increase/decay parameter for speciation. The control files containing important information for BAMM such as prior parameters and chain length were set individually for each tree using parameters estimated using the BAMMtools package (Rabosky et al. 2014), and all BAMM analysis were run in parallel using the GNU parallel shell tool (Tange 2011) and the latest BAMM version (2.5; Rabosky 2014). After running BAMM, the resulting files were imported in R and processed to retrieve the rates using the package BAMMtools. We note that the scenario with only varying speciation (SPvar) yields the most direct comparisons between the methods, given that the current version of BAMM does not allow a test for varying extinction rates. Even though the varying extinction scenario (EXvar) is only implemented in RPANDA, we decided to also use BAMM to estimate rates under this scenario to compare the methods, but more importantly to understand the behavior of both methods in a potentially more complicated scenario. This is justified given that it is an evolutionary plausible scenario, and from a given molecular phylogeny there is no a priori justification to exclude such scenario as a potential one. Empirical trees — We fitted both RPANDA and BAMM to 214 empirical trees. The trees used here are the same as the trees used by Lewitus and Morlon (2016), and include trees of all five major vertebrate groups (fishes, amphibians, reptiles, birds, and mammals). For RPANDA, we fit all four models used in the simulation analysis plus a fifth model that represented an exponential increase in both speciation and extinction (hereafter called BOTHexp). This model was chosen because it includes an additional scenario of expanding diversity (which might be plausible in empirical phylogenies). Depending on the parameter estimates it can represent either quasiconstant rates (when the increase parameter is close to 0) or increasing rates (when the increase parameter is different from 0), therefore providing a versatile scenario for testing. For BAMM, we estimated the set of parameters for each tree individually using the BAMMtools package for R. Since there was only one large tree (Muridae, 680 species) for which the “expected number of shifts” parameter should be set different than 1, we used 1 for all trees. All the estimated rates were then summarized following the same protocol as for the simulated trees, the only addition being that for the empirical trees, when using RPANDA, we estimated the values for each parameter via model averaging, in addition to analyzing the results for the best selected model. These estimates consisted in the parameter values in each of the four models averaged by their respective Akaike weights. We do not know which model is the true model that shaped the diversity of each empirical clade, whereas for the simulated data set we knew the true model in each scenario. Therefore, we used the model averaging approach only for the empirical estimates. Lastly, we estimated the correlation between the estimated rates of both methods (speciation, extinction, and net diversification) both at the tips and at the root using linear models to help us better understand the differences between the parameter estimates. IV: Performance analysis — All the parameter estimates for all 16,000 trees were tabulated and used to evaluate the performance of the two methods using Pearson’s correlation tests (Table 5). All parameter estimates shown in main results correspond to the best model for RPANDA, regardless of whether or not it is the single best model. Also, to explore whether recent studies using these methods have uncovered clades in decline, we surveyed the recent literature by analyzing the results from papers that cited either Morlon et al. (2011) and/or Morlon et al. (2016), and Rabosky (2014). Results Simulations Our results indicate that the collection of viable trees allowed us to greatly expand the parameter space explored by Morlon et al. (2011) and explore the total parameter space delimited by the values present in Table 1 in a comprehensive manner (Fig. 2a, b). Our simulations also resulted in trees with similar size and depth distributions for both scenarios (Fig. 2c, d), with trees ranging from 2 to more than 6000 species, and from 3.3 to 100 million years of age. Hence, the simulations of both decline scenarios resulted in a virtually indistinguishable distribution of phylogenies with respect to total time and richness. The branching patterns on the other hand are very different. The standardized gamma statistic distributions of both scenarios have very little overlap (Fig. 2e), being predominantly negative for SPvar and positive for more than 98% of the trees at EXvar. Fitting We present the results of three out of the four time slices (80rise, 80left, and 20left, respectively). The results for the 50left time slice is only presented in the Supplementary Material available on Dryad at http://dx.doi.org/10.5061/dryad.tt10r84 (Supplementary Figs. S1–S4 available on Dryad) because they are qualitatively similar, and intermediate to the 80left (discrete decline) and 20left (pronounced decline). The predominantly positive estimated values of net diversification rates for the 80rise time slice indicate that there is no bias towards falsely detecting decline for both methods in both scenarios (Panels a and d of Supplementary Figs. S1 and S2 available on Dryad), nor any bias in estimating individual rates of speciation and extinction (Figs. 3 and 4a, d, g, and j) during the diversity expansion phase. Hence we focus on the two other time slices, which represent our main interest in diversity decline. Figure 3. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the SPvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). Figure 3. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the SPvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). Figure 4. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the EXvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). Figure 4. View largeDownload slide Estimated versus simulated speciation rates at present (panels a–f) and estimated versus simulated extinction rates at present (panels g–l) for the EXvar scenario for RPANDA (panels a, b, c, g, h, and i) and BAMM (panels d, e, f, j, k, and l). Estimates are for the 80rise time slice (left column), 80left time slice (middle column), and 20left time slice (right column). For all RPANDA panels (a, b, c, g, h, and i), colors denote the best model selected, and in all panels the solid line denotes the ‘perfect fit’ (estimated rate equal to the simulated rate). The |$y$|-axes were rescaled to values that contained at least 95% of the trees for better visualization. (Colored figures in the online version; the summary of model selection by RPANDA can be found in table 2). SPvar — RPANDA selects the “true” model for the vast majority of the trees (91% for the 80left time slice and 87.5% for the 20left time slice—Table 2). The good performance of RPANDA on estimating net diversification is due to the nonbiased and accurate estimation of both speciation and extinction rates at present (Fig. 3b, c, g, h, and i and Supplementary Figs. S1, S3, S4, and S5 and Table S5 available on Dryad), with a very small bias for both rates in the smaller trees (Supplementary Figs. S3 and S4 available on Dryad). Extinction rate is frequently estimated to be very close to zero when the BOTHcst model is selected. Thus, RPANDA fails to detect decline for a few trees in both time slices, especially when the best model selected is BOTHcst. RPANDA estimates of net diversification show little bias for both time slices (80left and 20left—Supplementary Figs. S1 and S5, and Table S5 available on Dryad). When the most complex model (BOTHvar) was selected, the estimates of net diversification were highly negative, overestimating the pace of diversity decline for those trees (Supplementary Fig. 1b, c available on Dryad), especially for moderately sized trees. Table 2. Number of trees for which RPANDA choses each of the four tested models as the best model in both scenarios. Top level headings indicate the true simulated model in each scenarios. Numbers in bold show the instance where the correct model was chosen SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) Table 2. Number of trees for which RPANDA choses each of the four tested models as the best model in both scenarios. Top level headings indicate the true simulated model in each scenarios. Numbers in bold show the instance where the correct model was chosen SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) SPvar EXvar Model 80rise 80left 20left 80rise 80left 20left BOTHcst 335 (16.75%) 79 (3.95%) 108 (5.4%) 673 (33.65%) 222 (11.1%) 443 (22.15%) EXvar 69 (3.45%) 1 (0.05%) 3 (0.15%) 10 (0.5%) 10 (0.5%) 34 (1.7%) SPvar 1538 (76.9%) 1821 (91.05%) 1749 (87.45%) 1210 (60.5%) 1635 (81.75%) 1376 (68.8%) BOTHvar 58 (2.9%) 99 (4.95%) 140 (7%) 107 (5.35%) 135 (6.75%) 147 (7.35%) In the SPvar scenario, RPANDA detected a single best model for the majority of trees (Table 3). Additionally, it is worth noticing that even when more than one model were equally likely, the true model was among the best models for the vast majority of trees (Supplementary Table S1 available on Dryad). Similarly, BAMM is able to detect diversity decline for most trees in both time slices (Supplementary Fig. S1e, f available on Dryad), as well as a drop in speciation for at least 96.65% of trees (see Table 4). However, BAMM produced slightly negatively biased estimates for net diversification (Supplementary Fig. S1e, f and S8 available on Dryad). This bias is more evident for small trees (up to 100 tips), and also at the 80left time slice (Supplementary Fig. 8 available on Dryad). Interestingly, both speciation and extinction rates seem to be overestimated for trees up to 1000 tips, and this bias is seen more clearly at the 20left time slice although the number of trees is small and hence biases harder to infer (Fig. 3e, f, k, and l, and Supplementary Figs. S6 and S7 available on Dryad). Surprisingly, our results indicate that this bias seems to be stronger for speciation rates than for extinction rates (Supplementary Figs. S6 and S7 available on Dryad). Nevertheless, the biases in both rates seem to be coupled, that is, the overestimation in both rates for the same tree is proportional, since the net diversification rates are much less biased than each rate alone (Supplementary Fig. S8 available on Dryad). Table 3. Number of trees for which RPANDA choses each number of equally likely models based on the AICc differences (⁠|$\Delta$|AICc|${}>{}$|2) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) Table 3. Number of trees for which RPANDA choses each number of equally likely models based on the AICc differences (⁠|$\Delta$|AICc|${}>{}$|2) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) SPvar EXvar Number of equally likely models 80rise 80left 20left 80rise 80left 20left 1 1282 (64.1%) 1466 (73.3%) 1601 (80.05%) 856 (42.8%) 1125 (56.25%) 1274 (63.7%) 2 641 (32.05%) 530 (26.5%) 368 (18.4%) 912 (45.6%) 795 (39.75%) 637 (31.85%) 3 73 (3.65%) 3 (0.015%) 18 (0.9%) 190 (9.5%) 68 (3.4%) 77 (3.85%) 4 4 (0.02%) 1 (0.05%) 13 (0.65%) 42 (2.1%) 12 (0.6%) 12 (0.6%) Table 4. Number of trees for which BAMM infers decrease in speciation rates in each simulated scenario 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) Table 4. Number of trees for which BAMM infers decrease in speciation rates in each simulated scenario 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) 80rise 80left 20left SPvar 76 (96.2%) 67 (96.65%) 13 (99.35%) EXvar 246 (87.7%) 342 (82.9%) 509 (74.55%) The amount of extinction experienced by the clade had a differential effect on methods. For example, the comparison between the 80left and the 20left time slices suggests that RPANDA estimates were largely robust to increased amounts of extinction. The method showed a similar ability to detect negative diversification in both time slices (although the number of trees for which the method chooses the simplest model—BOTHcst—increases at the 20left time slice—Table 2). Conversely, BAMM estimates seemed to benefit from an increased duration of the diversity decline phase and detected diversity decline more frequently at 20left than at 80left time slice (Supplementary Fig. S1e, f available on Dryad). For both methods, the estimated speciation and extinction rates at the root of the trees (reconstructed values for RPANDA and estimated values for BAMM) are similar to the simulated values, regardless of the time slice for both methods (Supplementary Fig. S9 and Table S5 available on Dryad). RPANDA tends to overestimate speciation rates at the root more than BAMM, probably as a consequence of those rate estimates being influenced by the propagation of the error in the estimates of the rate at the present and the decay/increase parameter. This overestimation is greater when the best model selected by RPANDA has variation in extinction rates (BOTHvar). Accordingly, estimated extinction rates at the root are underestimated when the BOTHvar model is selected, being very close to 0. Table 5. Estimated |$R^{2}$| for the linear regression between rate estimates of BAMM and RPANDA for each rate and at each point in time (tip or root). Tip represents estimates at present and “root” at the start of the history of the clades. Values in bold are significant (⁠|$P<0.05$|⁠), and values in both bold and italic are considered poor associations regardless of being significant (⁠|$P<0.05$|⁠). Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 Table 5. Estimated |$R^{2}$| for the linear regression between rate estimates of BAMM and RPANDA for each rate and at each point in time (tip or root). Tip represents estimates at present and “root” at the start of the history of the clades. Values in bold are significant (⁠|$P<0.05$|⁠), and values in both bold and italic are considered poor associations regardless of being significant (⁠|$P<0.05$|⁠). Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 Lambda Mu Net diversification Root Tip Root Tip Root Tip Full data set |$-$|0.004 0.278 0.418 0.009 |$-$|0.003 0.23 Without BOTHexp |$-$|0.006 0.876 0.602 0.012 |$-$|0.004 0.33 Without varying extinction 0.629 0.861 0.692 0.7 0.433 0.239 Only varying extinction |$-$|0.009 0.05 0.019 0.018 0.003 0.122 EXvar — The performance of both methods was worse for the varying extinction scenario (EXvar). RPANDA rarely selected the true model as the best model (Table 2); at both the 80left and 20left time slices SPvar was the preferred model for most trees (Table 2), followed by the model with both rates constant. It is also worth noting that the method failed to detect the true simulated model for about 99.6% of the trees at the 80left time slice and for about 98.3% at the 20left time slice. Additionally, in the EXvar scenario RPANDA detects a single best model for much fewer trees than in the SPvar scenario (Table 3). Moreover, for the trees for which one or two models were selected as best (⁠|$\Delta$|AICc |$>$|2), the true (simulated) model was almost never among them (Supplementary Table S1 available on Dryad). Lastly, the estimates under the true model in this scenario (EXvar) were biased towards underestimation for both speciation and extinction at present (Supplementary Fig. S10g–l available on Dryad). This underestimation was more severe for extinction rate, and these biases combined made estimates of net diversification rates at the present to be positive under the true model for the majority of the trees (Supplementary Fig. S10i, l available on Dryad). Overall RPANDA seemed unable to detect varying extinction rates. Although RPANDA detects diversity decline for the EXvar scenario, it did by selecting the wrong model (see Table 2). It is interesting to note that RPANDA showed two distinct biases for net diversification rate estimates: it either estimated net diversification rates to be considerably more negative than the simulated values (regardless of tree size), or it tended to estimate this rate to be a lot less negative, or even positive, for a considerable amount of trees, especially those trees up to 1000 tips (Supplementary Fig. S13 available on Dryad). Nevertheless, the estimates for both speciation and extinction rates seemed nonbiased when we looked at all estimates together (Supplementary Figs. S11 and S12 available on Dryad). This implies that for any given tree RPANDA estimated at least one of the two rates at present to be considerably different from the simulated one. Like RPANDA, BAMM was able to detect negative net diversification rates for most trees, even when BAMM’s implemented model (time-varying speciation with fixed extinction) conflicted with the model used in simulation. However, BAMM estimates of net diversification rates are slightly less biased in all tree sizes (Supplementary Fig. S16 available on Dryad). For BAMM the bias is unidirectional, where estimates of net diversification rates tended to be less negative than the simulated values (Supplementary Fig. S16 available on Dryad). This bias difference is clearer for small trees. Speciation and extinction rates were slightly underestimated at the 20left time slice for trees up to 100 tips (Supplementary Figs. S14 and S15 available on Dryad). Furthermore, the rate of suggesting drops in speciation rates are higher than in the SPvar scenario (at least 74.55% of trees in the two scenarios—Table 4). Overall both methods are more likely to detect diversity decline at the 20left time slice than at the 80left time slice. Interestingly RPANDA detected negative diversification more often than BAMM at both time slices. At 80left time slice, RPANDA detected diversity decline for 1799 trees (89.95%), and BAMM detected decline for 1544 trees (77.2%). At the 20left time slice, RPANDA detected diversity decline for 1840 trees (92%) and BAMM for 1684 trees (84.2%) (Supplementary Fig. S2 available on Dryad). In the EXvar scenario, speciation and extinction rates at present are reasonably well estimated in both methods (represented by the small deviations from the perfect fit line in Fig. 4 and Supplementary Table S5 available on Dryad) although the net diversification rate deviates considerably from the simulated values. However, both speciation and extinction rates at the root do not match the simulated values (Supplementary Figs. S17 and S18 and Table S5 available on Dryad). In fact, both methods infer a considerable drop in speciation and a relatively constant extinction at both time slices for most trees, even though the simulated model had constant speciation rate and variable extinction. Even when RPANDA selected variable extinction as the generating model (145 trees for 80left and 181 trees for 20left—adding up trees for which EXvar and BOTHvar models were selected as best models in Table 2), it rarely provided good estimates of the extinction rate at the root of the tree (Supplementary Fig. 18 available on Dryad); on the other hand, BAMM does not allow for extinction to vary, which forces the method to estimate root extinction to have the same values than in the present. Surveying the literature with Google Scholar we found that, until January 2018, 149 empirical studies cited Rabosky (2014) and applied BAMM, whereas 24 empirical studies cited Morlon et al. (2011) and/or Morlon et al. (2016) and applied RPANDA. Several of these studies, however, provide only speciation rates or shift plots without rate values, rendering impossible to know if diversity declines have been detected. From those studies providing enough information for us to know whether or not they detected decline, 9.4% (14 studies) and 20.8% (5 studies) detected decline using BAMM or RPANDA, respectively (for the complete list of publications, see Supplementary Table S1 available on Dryad). The tabulation using RPANDA, although more favorable in terms of finding diversity decline, represent a smaller sample size. Analysis of Empirical Trees RPANDA indicates a decline in diversity (net diversification smaller than |$-$|0.01) for 31% of the analyzed clades, whereas BAMM detects decline for only 6% of clades from the same data set (see also Supplementary Fig. S19 available on Dryad). For only 3.7% of the trees, a negative net diversification rate was inferred for both methods (see also Supplementary Fig. S19 available on Dryad). Our results also indicate that most empirical trees for which decline was inferred also show negative corrected gamma values (Supplementary Fig. S24 available on Dryad). Speciation rates at present were estimated to be similar between both methods (Fig. 5a and Table 5), especially after removing the trees that RPANDA fitted the BOTHexp model, which produced very high, and sometimes nrealistic, extinction rate estimates (Table 4). One should note that the surplus of trees with negative net diversification rates when using RPANDA comes, from the most part, from those trees that fitted models with increase in both rates, the BOTHexp model (43 out of 67 trees; see also Supplementary Table S2 available on Dryad). Extinction rate estimates at present are not significantly correlated between methods, even after excluding the model BOTHexp (see Fig. 5, see “without BOTHexp” results for “Tip” on Table 5). Although RPANDA suggested a broader variation in extinction rates than BAMM (Fig. 5), many of those rates were estimated to be very small at present (e.g., 113 trees with extinction estimates |$\leqslant$|0.01). BAMM, on the other hand, very rarely estimated extinction rates to be so close to 0 (e.g., 2 trees with extinction estimates |$\leqslant$|0.01). The correlation between the estimates of net diversification at present are, as expected, somewhere in between those of the individual rates (Supplementary Fig. S19 available on Dryad; Table 5). Figure 5. View largeDownload slide Estimated rates for BAMM (⁠|$y$|-axis) and RPANDA (⁠|$x$|-axis) for (a, d) speciation rates at the tips and the root, respectively, and for (b, c, e, f) extinction rates at the tips and the root, respectively. The |$x$|-axis of panels b and e were rescaled from panels c and f, respectively, for better visualization, since there are estimated extinction rates from RPANDA that are greater than 10 events/lineage*myr (up to more than 1500 in the detail of c). The red line represents the perfect correlation (estimates from BAMM and RPANDA are equal), and the blue line represents a linear fit between the two variables. |$R^{2}$| values are shown in Table 5. Speciation rates are similar between the two methods, whereas extinction can be quite distinct between the two methods. (Colored figures in the online version). Figure 5. View largeDownload slide Estimated rates for BAMM (⁠|$y$|-axis) and RPANDA (⁠|$x$|-axis) for (a, d) speciation rates at the tips and the root, respectively, and for (b, c, e, f) extinction rates at the tips and the root, respectively. The |$x$|-axis of panels b and e were rescaled from panels c and f, respectively, for better visualization, since there are estimated extinction rates from RPANDA that are greater than 10 events/lineage*myr (up to more than 1500 in the detail of c). The red line represents the perfect correlation (estimates from BAMM and RPANDA are equal), and the blue line represents a linear fit between the two variables. |$R^{2}$| values are shown in Table 5. Speciation rates are similar between the two methods, whereas extinction can be quite distinct between the two methods. (Colored figures in the online version). Inference of speciation rate dynamics showed a remarkable difference between the two methods. RPANDA suggests that the vast majority of empirical trees (166 out of 214—77%) fitted models where speciation was either constant (123 out of 214—57%) or increasing (43 out of 214—20%) (Fig. 6a). BAMM, on the other hand, suggested that the majority of empirical trees (160 out of 214—75%) showed a relative decrease in speciation (Fig. 6a; most points below zero on the |$y$|-axis). Surprisingly, estimates of extinction rate at the root are more correlated than at the present (Fig. 5e and Table 5). This is partially due to the restriction of initial values related to the particularities of two models (SPvar, BOTHcst) that are implemented in BAMM, and to the fact that three (EXvar, BOTHvar, and BOTHexp) out of the five tested models in RPANDA, are strongly constrained to start with very low (or zero) extinction rates (see Fig. 5). Figure 6. View largeDownload slide Relative speciation rate variation, calculated as [lambda at tips |$-$| lambda at root]/max(lambda at tips, lambda at root) of RPANDA as a function of the same index for BAMM. Positive values indicate an increase in speciation towards the present, whereas negative values indicate decreasing speciation rates through time; 0 indicates constant rates. a) Points are colored according to which model detected decline for each empirical tree; b) Points are colored according to which model was selected as the best model by RPANDA. It is possible to notice that the points that deviate the most from the perfect association (solid line) have models with varying extinction rates as the best model, and that the points that represent a model that is the same as the model used by BAMM fall very close to the red line (green dots on the bottom left part of panel b). (Colored figures in the online version). Figure 6. View largeDownload slide Relative speciation rate variation, calculated as [lambda at tips |$-$| lambda at root]/max(lambda at tips, lambda at root) of RPANDA as a function of the same index for BAMM. Positive values indicate an increase in speciation towards the present, whereas negative values indicate decreasing speciation rates through time; 0 indicates constant rates. a) Points are colored according to which model detected decline for each empirical tree; b) Points are colored according to which model was selected as the best model by RPANDA. It is possible to notice that the points that deviate the most from the perfect association (solid line) have models with varying extinction rates as the best model, and that the points that represent a model that is the same as the model used by BAMM fall very close to the red line (green dots on the bottom left part of panel b). (Colored figures in the online version). Discussion Although the development of new methods (Morlon et al. 2011; Rabosky 2014) might in principle allow estimates of net diversification rates to be negative, the question of whether extinction dynamics can be properly estimated from molecular phylogenies remains contentious (Nee et al. 1994; Purvis 2008; Quental and Marshall 2010; Rabosky 2010; Morlon et al. 2011; Beaulieu and O’Meara 2015; Rabosky 2016; Moore et al. 2016; Rabosky et al. 2017). Our results shed new light on the performance of these methods by considering a wide range of decline scenarios. When looking at simple simulated scenarios, we find that both methods are able to detect decline at the present, but both also fail to characterize the temporal dynamics when extinction rate varies. On the other hand, when looking at empirical phylogenies, we find that estimates from both methods are in general not correlated. By the virtue of knowing the true dynamics, simulation studies are extremely valuable to investigate our ability to detect negative diversification rates based solely on molecular phylogenies, but to date this has been limitedly done to only one of the methods discussed here (Morlon et al. 2011). Simulation Perspective on Our Ability to Infer Negative Diversification Rates Although Morlon et al. (2011) suggested that RPANDA was able to detect diversity decline, their study was limited in several important respects: 1) their parameter space was restrictive when compared to broader ranges of biologically plausible values of speciation and extinction rates; 2) it did not explicitly consider method performance at different intensities of decline (here represented by the different time slices); 3) it only examined performance of a single method. Here, we not only explicitly test the ability of BAMM to detect negative diversification rate for the first time, but we expanded the findings by Morlon et al. (2011) by exploring broader parameter space and different points in time of a given diversification scenario. Unexpectedly, we now show that even when the underlying rate dynamics is wrong (wrong model in the RPANDA analysis; or wrong overall dynamics in BAMM), both methods are able to infer negative net diversification at present, at least for the scenarios simulated here. Moreover, our results indicate that the different phase of a tree’s history (our scenarios at 80left, 50left, and 20left time slices) has only a slight impact on the ability of both models in detecting diversity decline at present in the both simulated scenarios (Supplementary Table S5 available on Dryad). This is surprising given that the phylogenetic overall signature of diversification dynamics (the branching patterns) changes as the diversification processes unfolds and the clade ages, a pattern previously revealed by changes in the gamma statistic (Liow et al. 2010; Quental and Marshall 2011). It should be noted that our results also suggested that inferring reliable negative net diversification rates at the recent past is not a proxy for the reliability of the whole diversification dynamics reconstruction. Can We Properly Infer the Diversification Dynamics of Clades in Decline? Our results suggest that the accuracy of the methods in depicting the temporal dynamics, either by choosing the best model (in the case of RPANDA) or by estimating the rates at both the present and at the root of the trees (in the case of BAMM) varied drastically according to the scenario explored. We clearly show that when extinction was allowed to vary (EXvar scenario), both methods performed worse than when extinction was held constant (SPvar scenario) (Supplementary Table S5 available on Dryad). Given that RPANDA is not limited by preimplemented models (e.g., in our analysis we included varying extinction rate models), its inability to infer varying extinction rates suggests that this poorer performance of both methods might be more related to the limited information inherent to molecular data (e.g., no direct information on extinction) than to the method or implementation used. Under the EXvar scenario, speciation rate estimates at the root (Supplementary Fig. S17 available on Dryad) are significantly higher than the true values for both BAMM and RPANDA. Similarly, both methods suggested a significant drop in speciation rate (Fig. 4 and Supplementary Fig. S17 available on Dryad) when in reality speciation was held constant. Moreover, the prevalence of models chosen by RPANDA with varying speciation in the simulated EXvar scenario indicates that RPANDA detected variation in net diversification but explained it as dynamically changing speciation rate. Given that those methods rely on indirect inference of the past, the amount of information available for both methods becomes thinner as we move towards the root of the tree. In fact, both RPANDA and BAMM showed greater associated error for parameters at the root regardless of the simulated scenario (Supplementary Figs. S9 and S18 available on Dryad), a pattern similar to what was found by Moore et al. (2016). This might not be surprising but we argue that it represents one important aspect of the two methods and (likely) the nature of signal in phylogenetic data sets. Methods based on molecular phylogenies could be seen as binoculars with limited range: they allow researchers to infer important aspects of the recent history of a clade, but it loses power and focus as one tries to look towards the deep past. It is worth noting, the gamma statistic (Pybus and Harvey 2000) for trees under the EXvar scenario are virtually all positive (Fig. 2), as expected by previous simulation studies (Pybus and Harvey 2000; Rabosky and Lovette 2008; Quental and Marshall 2009). Although gamma itself is not informative about whether a clade is in decline, our simulations suggest that the gamma statistic when combined with RPANDA and BAMM results could be helpful in discerning among decline scenarios. For instance, if BAMM and RPANDA suggest a decrease in speciation rates through time, but the gamma value is positive, this might in fact represent an increase in extinction rates instead of the variation in speciation (see Results section, Fig. 2). Comparing the Performance of Both Methods in Empirical Phylogenies Although both methods performed similarly (either well or badly) in our simulation scenarios, it is interesting to note that both inferred very different diversification dynamics when applied to our collection of 214 empirical trees. It is especially interesting to note that both methods only agree with respect to negative diversification rates on about 3.7% of the trees. Of course we do not know the true underlying dynamics for the empirical trees, but the comparison is informative nonetheless because a lack of convergence signals that one (or both) method(s) is providing incorrect estimates. Our empirical analysis also reinforces the idea that inferences at deep time might be more challenging than inferences at the recent present. Both methods converged in similar speciation rate estimates at present but inferred very different speciation rates at the root (Fig. 5). BAMM tends to suggest an overwhelming drop in speciation towards the present, irrespective if the lineages were in diversity decline or not while RPANDA show many different patterns of rate variation. This discrepancy might be driven by the important difference on how flexible the two methods are regarding modeling extinction rates because most discrepant estimates in speciation dynamics come from those cases where RPANDA inferred changes in extinction (Fig. 5b, c). Extinction rate estimates at the presents showed an even stronger discrepancy across the methods (Fig. 5b, c). This is not surprising given that extinction rates can never vary with time in BAMM, but this discrepancy has important consequences when it comes to detecting diversity decline. Depending on the method used, one would have a very different view of how common is diversity decline in the real world (31% vs 6% of clades in diversity decline, for RPANDA and BAMM, respectively). This casts serious doubts on our current ability to estimate extinction rates (Rabosky 2010) and, more importantly, detect diversity decline based on current available methods. Our analysis of empirical data sets indicates a conflict between the speciation dynamics for the set of analyzed trees and results from some prior studies. In our study, results from RPANDA suggest that increasing speciation rates are rather common. On the other hand, both the gamma statistic (McPeek 2008; Phillimore and Price 2008) and a more complex model fitting study (Morlon et al. 2010) suggest that slowdowns are pervasive. It should be noted that the method used in Morlon et al. (2010) did not allow (at least in that particular paper) extinction rates to be higher than speciation rates. RPANDA (Morlon et al. 2011) allows net diversification to be negative, hence our results for the empirical data set lead us to infer that either a drop in speciation might be a common phenomenon and RPANDA misfit speciation dynamics for clades in diversity decline; or that a drop in speciation is not so prevalent as once thought. Lastly, the association between underlying drivers of diversity decline estimated by BAMM and RPANDA and the gamma values might be worth exploring. For BAMM the vast majority of clades in decline show negative gamma values (which would suggest a decrease in speciation rates), and those show a drop in speciation. This association between negative gamma values and a drop in speciation is not as strong for the RPANDA estimates. Whereas most clades in decline show negative gamma values, RPANDA did not assigned scenarios of decrease in speciation for many of them (see empirical summary in Supplementary Material available on Dryad). Those that show inconsistency between gamma values and the role of speciation should be seen with caution (Quental and Marshall 2009, 2011). The literature survey revealed a lack of empirical examples of negative diversification rates in the literature that is quite surprising. Even though this paucity could be related to a lack of studies using such methods, the total number of papers analyzed in our literature survey suggests otherwise. Another potential explanation for this empirical lack of negative diversification rates could be that empirical molecular phylogenies are either very recent radiations, or much more heterogeneous than the ones we typically simulate. If this were the case, the prevalence of molecular phylogenies not showing diversity decline would not be surprising, but rather expected. Empirical trees could for example include a radiating clade that masks a potential diversity decline of a subclade (see Moyne and Neige 2007; or the contrast between Morlon et al. 2011 and Rabosky 2014 for potential examples). We cannot evaluate this for our literature survey, but for the 214 empirical trees analyzed here this does not seem to be the case. Our BAMM analysis did not find any significant shifts for the vast majority of those empirical trees (see Results section). Additionally the empirical trees analyzed here have crown group ages that range from 6.7 myr to 131 myr, with a median value of about 30 myr. It is possible that most of the recent radiations are indeed not in decline, but judging from what we know from the fossil record (e.g., Silvestro et al. 2015; Pires et al. 2017), we would expect many of those older radiations to have had enough time to pass their expanding phase. Although we cannot use the empirical results to evaluate methods performance, their lack of convergence casts serious doubts on our ability to properly estimate the underlying diversification dynamics, especially at deep time. Final Remarks: How Reliable are Estimates of Diversification Dynamics The last decades have witnessed the constant development of new comparative methods (Morlon et al. 2011; FitzJohn 2012; Rabosky 2014), and the turmoil of optimistic and pessimistic views on how well can we infer diversification dynamics solely from molecular phylogenies. This field has recently arrived to its final frontier, the idea that one could detect diversity decline, a scenario where the recent is experiencing more extinction than speciation even though no direct information on extinction is stored in phylogenies. Only one study (Morlon et al. 2011) has explicitly tested the ability of phylogeny-based models on detecting negative diversification rates, and only a few empirical studies have in fact detect this scenario. For the first time to our knowledge, we explicitly tested the ability of BAMM (Rabosky 2014) to detect negative diversification rates. Although our simulation results suggest that these methods might be able to detect negative diversification rates, we also showed that the inferred dynamics (either using BAMM or RPANDA) might be utterly wrong. Additionally, our simulation results and the lack of convergence in the empirical analysis suggested that the empirical pattern of drop in speciation previously thought to be prevalent (Morlon et al. 2010) might in fact result, at least in part, from model “misbehavior.” Hence, we suggest that current methods might have a difficult time on simultaneously inferring deep time dynamics and recent dynamics, especially for clades experiencing diversity decline. Furthermore, based on our results we reinforce that researchers should try to use external sources of information (e.g., fossil record when available) when trying to reconstruct diversification dynamics in deep time. This might be especially relevant if there is external evidence for considerable variation in extinction rates. Although we remain pessimistic about our ability to properly infer the diversification dynamics in deep time, especially for clades experiencing diversity decline, we believe that our results might help map the scenarios where those models are more likely to succeed. Lastly, and perhaps more optimistically, we suggest that current phylogenetic methods might at least properly infer very recent dynamics. Supplementary material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.tt10r84. Funding This work was supported by São Paulo Research Foundation (FAPESP) [grants #2014/03621-9 to G.B., #2012/02038-2 to L.R.V.A. and #2012/04072-3 to T.B.Q.]; and NSF [DEB-0918748 to M.E.A.]. Acknowledgments We would like to thank Dr. Hélène Morlon and Dr. Daniel Rabosky for their helpful comments with RPANDA and BAMM. We further thank Dr. Hélène Morlon for all the assistances throughout the whole development of the article. We also thank Mericien Venzon for the aid in running BAMM. We extend our gratitude to Dr. Eric Lewitus and Dr. Hélène Morlon for kindly providing the empirical trees and all students at LabMeMe (Laboratory of Macroevolution and Macroecology) at the University of São Paulo, Dr. Luke Harmon, and two anonymous reviewers for their valuable comments on the study. References Alencar L.R. , Quental T.B. , Grazziotin F.G. , Alfaro M.L. , Martins M. , Venzon M. , Zaher H. 2016 . Diversification in vipers: phylogenetic relationships, time of divergence and shifts in speciation rates . Mol. Phylogenet. Evol. 105 : 50 – 62 . Google Scholar Crossref Search ADS PubMed Alroy J. 1996 . Constant extinction, constrained diversification, and uncoordinated stasis in North American mammals. Palaeogeogr. Palaeoclimatol. Palaeoecol . 127 ( 1–4 ): 285 – 311 Google Scholar Crossref Search ADS Alroy J. 2010 . Fair sampling of taxonomic richness and unbiased estimation of origination and extinction rates. In: Alroy J., Hunt G., editors. Quantitative methods in paleobiology . The Paleontological Society Papers. Vol. 16. The Paleontological Society . Bethesda, MD, USA. 316pp . Bapst D.W. 2012 . paleotree: an R package for paleontological and phylogenetic analyses of evolution . Methods Ecol. Evol. 3 5 : 803 – 807 . Google Scholar Crossref Search ADS Barraclough T.G. , Nee S. 2001 . Phylogenetics and speciation . Trends Ecol. Evol. 16 7 : 391 – 399 . Google Scholar Crossref Search ADS PubMed Beaulieu J.M. , O’Meara B.C. 2015 . Extinction can be estimated from moderately sized molecular phylogenies . Evolution 69 : 1036 – 1043 . Google Scholar Crossref Search ADS PubMed Burin G., Kissling W.D., Guimarães P.R. , Şekercioğlu Ç.H. , Quental T.B. 2016 . Omnivory in birds is a macroevolutionary sink . Nat. Commun. 7 : 11250 . Google Scholar Crossref Search ADS PubMed Burnham K.P. , Anderson D. 2002 . Model selection and multi-model inference. A practical information-theoretic approach. New York, NY, USA : Springer . ISBN 978-0-387-22456-5. FitzJohn R.G. 2012 . Diversitree: comparative phylogenetic analyses of diversification in R . Methods Ecol. Evol. 3 : 1084 – 1092 . Google Scholar Crossref Search ADS Foote M. 2007 . Symmetric waxing and waning of marine invertebrate genera . Paleobiology 33 4 : 517 – 529 . Google Scholar Crossref Search ADS Gilinsky N.L. Bambach R.K. 1987 . Asymmetrical patterns of origination and extinction in higher taxa . Paleobiology 13 : 427 – 445 . Google Scholar Crossref Search ADS Heath T.A. , Huelsenbeck J.P. , Stadler T. , 2014 . The fossilized birth-death process for coherent calibration of divergent-time estimated . Proc. Natl. Acad. Sci. USA 11 29 : E2957 – 66 . Google Scholar Crossref Search ADS Huang H. , Rabosky D.L. 2014 . Sexual selection and diversification: reexamining the correlation between dichromatism and speciation rate in birds . Am. Nat. 184 : E101 – E114 . Google Scholar Crossref Search ADS PubMed Lewitus E. , Morlon H. , 2016 . Natural constraints to species diversification . PLoS Biol. 14 8 : e1002532 . Liow L.H. , Stenseth N.C. , 2007 . The rise and fall of species: implications for macroevolutionary and macroecological studies . Proc. Biol. Sci. 274 1626 : 2745 – 2752 . Google Scholar Crossref Search ADS PubMed Liow L.H. , Quental T.B. , Marshall C.R. 2010 . When can decreasing diversification rates be detected with molecular phylogenies and the fossil record? Syst . Biol. 59 6 : 646 – 659 . Liow L.H. , Reitan T. , Harnik P.G. 2015 . Ecological interactions on macroevolutionary time scales: clams and brachiopods are more than ships that pass in the night. Ecol. Lett. 18 10 : 1030 – 1039 . Google Scholar Crossref Search ADS PubMed McPeek M.A. 2008 . The ecological dynamics of clade diversification and community assembly . Am. Nat. 172 6 : E270 – 84 . Google Scholar Crossref Search ADS PubMed Moore B.R. , Höhna S. , May M.R. , Rannala B. , Huelsenbeck J.P. 2016 . Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures . Proc. Natl. Acad. Sci. USA 113 34 : 9569 – 9574 . Google Scholar Crossref Search ADS Morlon H. 2014 . Phylogenetic approaches for studying diversification . Ecol. Lett. 17 : 508 – 525 . Google Scholar Crossref Search ADS PubMed Morlon H. , Parsons T.L. , Plotkin J.B. 2011 . Reconciling molecular phylogenies with the fossil record . Proc. Natl. Acad. Sci. USA 108 : 16327 – 16332 . Google Scholar Crossref Search ADS Morlon H. , Potts M.D. , Plotkin J.B. 2010 . Inferring the dynamics of diversification: a coalescent approach. PLoS Biol. 8 9 : e1000493 . Morlon H. , Lewitus E. , Condamine F.L. , Manceau M. , Clavel J. , Drury J. 2016 . RPANDA: an R package for macroevolutionary analyses on phylogenetic trees . Methods Ecol. Evol. 7 : 589 – 597 . Google Scholar Crossref Search ADS Moyne S. , Neige P. 2007 . The space-time relationship of taxonomic diversity and morphological disparity in the Middle Jurassic ammonite radiation . Palaeogeogr. Palaeoclimatol. Palaeoecol. 248 : 82 – 95 . Google Scholar Crossref Search ADS Nee S. , May R.M. , Harvey P.H. 1994 . The reconstructed evolutionary process . Philos. Trans. R. Soc. Lond. B Biol. Sci. 344 : 305 – 311 . Google Scholar Crossref Search ADS PubMed Nichols J.D. Pollock K.H. 1983 . Estimating taxonomic diversity, extinction rates, and speciation rates from fossil data using capture-recapture models . Paleobiology 9 2 : 150 – 163 . Google Scholar Crossref Search ADS Peters S.E. 2005 . Geologic constraints on the macroevolutionary history of marine animals. Proc. Natl. Acad. Sci. USA 102 35 : 12326 – 12331 . Google Scholar Crossref Search ADS Phillimore A.B. , Price T.D. 2008 . Density-dependent cladogenesis in birds. PLoS Biol. 6 3 : e71 . Google Scholar Crossref Search ADS PubMed Pires M.M. , Silvestro D. , Quental T.B. 2017 . Interactions within and between clades shaped the diversification of terrestrial carnivores . Evolution 71 7 : 1855 – 1864 . Google Scholar Crossref Search ADS PubMed Purvis A. 2008 . Phylogenetic approaches to the study of extinction . Annu. Rev. Ecol. Evol. Syst. 39 : 301 – 319 . Google Scholar Crossref Search ADS Pybus O.G. , Harvey P.H. 2000 . Testing macro-evolutionary models using incomplete molecular phylogenies . Proc. Biol. Sci. 267 : 2267 – 2272 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2009 . Extinction during evolutionary radiations: reconciling the fossil record with molecular phylogenies . Evolution 63 12 : 3158 – 3167 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2010 . Diversity dynamics: molecular phylogenies need the fossil record . Trends Ecol. Evol. (Amst.) 25 8 : 434 – 441 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2011 . The molecular phylogenetic signature of clades in decline . PLoS One 6 : e25780 . Google Scholar Crossref Search ADS PubMed Quental T.B. , Marshall C.R. 2013 . How the red queen drives terrestrial mammals to extinction . Science 341 : 290 – 292 . Google Scholar Crossref Search ADS PubMed R Core Team 2016 . R: a language and environment for statistical computing. Vienna, Austria : R Foundation for Statistical Computing. Available from: https://www.R-project.org. Rabosky D.L. 2010 . Extinction rates should not be estimated from molecular phylogenies . Evolution 64 : 1816 – 1824 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. 2014 . Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS One 9 2 : e89543 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. 2016 . Challenges in the estimation of extinction from molecular phylogenies: a response to Beaulieu and O’Meara . Evolution 70 1 : 218 – 228 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. , Lovette I.J. 2008 . Density-dependent diversification in North American wood warblers . Proc. Biol. Sci. 275 1649 : 2363 – 2371 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. , Grundler M. , Anderson C. , Title P. , Shi J.J. , Brown J.W. , Huang H. , Larson J.G. 2014 . BAMMtools: an R package for the analysis of evolutionary dynamics on phylogenetic trees . Methods Ecol. Evol. 5 7 : 701 – 707 . Google Scholar Crossref Search ADS Rabosky D.L. , Title P.O. , Huang H. 2015 . Minimal effects of latitude on present-day speciation rates in new world birds . Proc. Biol. Sci. 282 : 20142889 . Google Scholar Crossref Search ADS PubMed Rabosky D.L. , Mitchell J.S. , Chang J. 2017 . Is BAMM flawed? Theoretical and practical concerns in the analysis of multi-rate diversification models. Syst Biol. 66 : 477 – 498 . Google Scholar Crossref Search ADS PubMed Raup D.M. 1986 . Biological extinction in earth history . Science 231 : 528 – 1533 . Revolution Analytics , Weston S. , 2015a . foreach: ProvidesForeach Looping Construct for R. R package version 1.4.3. Available from: https://CRAN.R-project.org/package=foreach. Revolution Analytics, Weston S., 2015b . doMC: foreach parallel adaptor for ‘parallel’. R package version 1.3.4. Available from: https://CRAN.R-project.org/package=doMC. Ricklefs R.E. 2007 . Estimating diversification rates from phylogenetic information . Trends Ecol. Evol. 22 11 : 601 – 610 . Google Scholar Crossref Search ADS PubMed Roopnarine P.D. 2006 . Extinction cascades and catastrophe in ancient food webs . Paleobiology 32 1 : 1 – 19 . Google Scholar Crossref Search ADS Sakamoto M. , Benton M.J. , Venditti C. 2016 . Dinosaurs in decline tens of millions of years before their final extinction . Proc. Natl. Acad. Sci. USA 113 18 : 5036 – 5040 . Google Scholar Crossref Search ADS Sepkoski J.J. Jr. 1978 . A kinetic model of phanerozoic taxonomic diversity I . Analysis of marine orders. Paleobiology 4 3 : 223 – 251 . Google Scholar Crossref Search ADS Sepkoski J.J. Jr. 1979 . A kinetic model of phanerozoic taxonomic diversity II . Early phanerozoic families and multiple equilibria. Paleobiology 5 3 : 222 – 251 . Sepkoski J.J. Jr , 1984 . A kinetic model of phanerozoic taxonomic diversity . III. Post-paleozoic families and mass extinctions. Paleobiology 10 2 : 246 – 267 . Silvestro D. , Schnitzler J. , Liow L.H. , Antonelli A. , Salamin N. 2014 . Bayesian estimation of speciation and extinction from incomplete fossil occurrence data . Syst. Biol. 63 : 349 – 367 . Google Scholar Crossref Search ADS PubMed Silvestro D. , Antonelli A. , Salamin N. , Quental T.B. 2015 . The role of clade competition in the diversification of North American canids . Proc. Natl. Acad. Sci. USA 112 : 8684 – 8689 . Google Scholar Crossref Search ADS Slater G.J. , Harmon L.J. , Alfaro M.E. 2012 . Integrating fossils with molecular phylogenies improves inference of trait evolution . Evolution 66 : 3931 – 3944 . Google Scholar Crossref Search ADS PubMed Stadler T. 2013 . Recovering speciation and extinction dynamics based on phylogenies . J. Evol. Biol. 26 6 : 1203 – 1219 . Google Scholar Crossref Search ADS PubMed Stanley S.M. 2007 . An analysis of the history of marine animal diversity Paleobiology 33 ( sp6 ): 1 – 55 . Google Scholar Crossref Search ADS Starrfelt J. , Liow L.H. 2016 . How many dinosaur species were there? Fossil bias and true richness estimated using a Poisson sampling model. Philos. Trans. R. Soc. Lond. B Biol. Sci. 371 1691 : 20150219 . Tange O. 2011 . GNU parallel—the command-line power tool . login: USENIX Mag. 1 36 : 42 – 47 . Wickham H. , Francois R. 2016 . dplyr: a grammar of data manipulation. R package version 0.5.0 . Available from: https://CRAN.R-project.org/package=dplyr. © The Author(s) 2018. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. All rights reserved. For permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

### Journal

Systematic BiologyOxford University Press

Published: Jan 1, 2019

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations