The choice of tree prior and molecular clock does not substantially affect phylogenetic inferences of diversification rates

Introduction

Statistical comparative methods use phylogenetic trees to gain insight into macroevolutionary patterns and processes (Felsenstein, 1985; Harvey & Pagel, 1991; O’Meara, 2012; Rabosky, 2016; Harmon, 2018). Branch lengths and node ages provide information about the rate of lineage accumulation through time (Nee, May & Harvey, 1994; Nee, 2006; Ricklefs, 2007; Pyron & Burbrink, 2013) and are instrumental in characterizing the underlying processes generating global patterns of biodiversity (Schluter & Pennell, 2017). A typical workflow uses a point estimate of a phylogenetic tree or a distribution of trees to estimate macroevolutionary parameters, such as the rate of lineage accumulation (speciation) or extinction, which are often compared across groups to provide insight into diversification rates and the tempo of evolution (Nee, Mooers & Harvey, 1992; Magallón & Sanderson, 2001; Alfaro et al., 2009; Rabosky, 2014). However, parameter estimates are dependent on the tree from which they are inferred (Felsenstein, 1985). Most inference procedures assume that a tree is estimated without error, but, because branch lengths are fundamental to estimates of diversification parameters, uncertain phylogenies can be expected to yield uncertain estimates. Several studies have focused on the causes of parameter misestimation when fitting diversification models to trees (Nee, May & Harvey, 1994; Barraclough & Nee, 2001; Revell, Harmon & Glor, 2005; Cusimano & Renner, 2010; Rabosky, 2010; Moore et al., 2016), and a handful have evaluated uncertainty in phylogenetic estimation explicitly in the context of estimating diversification rates from phylogenies under specific simulation conditions (Revell, Harmon & Glor, 2005; Wertheim & Sanderson, 2011; Marin & Hedges, 2018).

Theoretical advances have expanded the scope of phylogenetic comparative methods for studying diversification. Historically, models only assumed a constant rate of lineage diversification or extinction (Nee, May & Harvey, 1994). More modern approaches utilize phylogenies to determine where and/or when shifts in the rates of speciation and extinction take place (see Pyron & Burbrink, 2013) or estimate rates that depend on species’ traits (Maddison, Midford & Otto, 2007; FitzJohn, Maddison & Otto, 2009; FitzJohn, 2010).

It has been shown that phylogenetic uncertainty and error in tree estimation can directly impact the results of diversification studies. For example, Revell, Harmon & Glor (2005) demonstrated that underparameterization of the model of nucleotide sequence evolution as part of the process of phylogenetic estimation can produce apparent slowdowns in the rate of diversification as quantified by Pybus & Harvey’s (2000) gamma statistic. Additionally, errors in branch lengths (Wertheim & Sanderson, 2011) and biased taxonomic sampling can both affect estimates (Höhna, 2014). Taken together, these studies suggest that phylogenetic error can affect the estimation of diversification rates.

Bayesian methods of phylogenetic inference produce posterior distributions of trees, and, therefore, diversification rates can be estimated across such distributions to quantify uncertainty. The use of Bayesian approaches in phylogenetics has increased in recent years due in part to the availability of software, including Bayesian Evolutionary Analysis by Sampling Trees (BEAST) (Drummond et al., 2012) and MrBayes (Ronquist et al., 2012). BEAST is a Java application that has seen widespread use in the phylogenetics community due to its ease-of-use, intuitive graphical user interface, and implementation of a number of phylogenetic and population genetic models. BEAST may also be run from the command line and can leverage GPU hardware, facilitating phylogeny reconstruction on high-performance computing architectures. Users can specify an analysis by passing options from the command line or through a GUI to a bundled application, BEAUti, which produces the XML input file required for BEAST. In a typical analysis, this XML defines models of sequence evolution, a choice of branching model (i.e., tree prior), and a choice of molecular clock, among other possible configurations.

The impact that the choice of priors governing the molecular clock and branching process (or “tree prior”) in molecular phylogenetics is an active area of research. Commonly used tree priors for inference among multiple species are the Yule (1925) and birth–death (BD; Kendall, 1948; Nee, May & Harvey, 1994; Gernhard, 2008; Stadler, 2013) models, whereas coalescent-based priors are suitable for phylogenetic and population genetic studies within a single species (Kingman, 1982; see Drummond et al., 2002, 2005). Here, we focus on the Yule and BD models. The Yule model is the simplest of a group of continuous-time branching processes; it has one parameter, λ, the instantaneous per-lineage rate of speciation, that is constant across the tree. The BD model is also a continuous-time process but includes a probability that a lineage will go extinct (and, therefore, leave no descendants). This model has two parameters, λ and μ, the instantaneous per-lineage rates of speciation and extinction, both of which are constant across the tree in their original characterizations. In practice, many approaches re-parameterize the model using r = (λ − μ) and ε = (μ/λ), the net diversification rate and relative extinction rate, respectively. In general, estimates of r have greater precision than ε (Nee et al., 1994; Nee, May & Harvey, 1994; FitzJohn, Maddison & Otto, 2009). Upon selecting these tree priors when using BEAST, a prior distribution (technically, hyperpriors in a hierarchical Bayesian context) must be specified on λ or on r and ε for Yule or BD, respectively.

Diversification rates can be estimated from phylogenies using likelihood-based approaches that rely on branching times (see Stadler, 2013). As a result, it is reasonable to assume that different branching models could have an impact on diversification rate estimates by virtue of altering branch lengths. Several studies have explored the impact of the tree prior on the resulting phylogenetic estimates. As part of an investigation of relaxed clock models, Ho et al. (2005) identified an impact of the choice of birth and death rate upper bounds in concert with the fraction of lineages sampled, particularly with respect to internal branches. Furthermore, in Ritchie, Lo & Ho (2017), the authors explore the impact of Yule and BD (and, additionally, coalescent) tree priors in the context of the multispecies coalescent to determine whether prior misspecification has an impact on phylogenetic accuracy. Through simulations and applications to empirical datasets, they concluded that phylogenies are not substantially affected by tree prior misspecification. However, node times may be influenced by the choice of prior in combination with among- and within-lineage sampling. Additionally, Brown & Yang (2010) found that, for shallow phylogenies, nodes depths are generally robust to the choice of prior. However, they concluded that a Dirichlet prior, in contrast to BD, produces more reasonable estimates as the depth of the phylogeny increases.

In addition to priors for branching process parameters, Bayesian phylogenetic analysis also requires the specification of a particular model for rates of evolution across the tree. For example, BEAST gives users the choice of using a strict (global) molecular clock or an uncorrelated log-normal (UCLN) relaxed molecular clock, among other flavors of local clocks (Drummond et al., 2012). The strict clock assumes a constant, global rate of sequence evolution across the tree (Zuckerkandl & Pauling, 1962), while the UCLN relaxed clock assumes branch-specific rates are drawn from a discretized log-normal distribution independently for every branch in the tree (Drummond et al., 2006). Hyperpriors are placed on the mean rate of evolution for the strict clock and the mean and standard deviation of the log-normal distribution for the UCLN relaxed clock. As the name implies, the UCLN molecular clock assumes that rates of evolution are not correlated among branches. This is in contrast with approaches commonly used to scale phylogenies after estimation, such as the penalized likelihood approach implemented in r8s (Sanderson, 2003) or treePL (Smith & O’Meara, 2012), which may inappropriately infer similar rates among closely-related lineages. However, the effect of selecting uncorrelated models over autocorrelated models may not always be clear and warrants further consideration (Ho et al., 2005; Lepage et al., 2007).

As with tree priors, the choice of molecular clock could also be expected to affect diversification rate estimates a priori by impacting branch lengths. Lepage et al. (2007) compare several clock models and show that clock choice can impact the estimates of divergence times. Furthermore, they find that clock misspecification can have a larger impact than the choice of branching prior. Previous work has also shown that relaxed clock models produce reasonable estimates of rates when substantial rate variation is not observed (Ho et al., 2005).

From the results outlined above, there is reason to believe that the choice of priors can affect the estimation of diversification parameters by virtue of altering the distribution of branch lengths. This has been explored specifically in several studies. For example, the effects of tree reconstruction on diversification rate estimates were studied by Wertheim & Sanderson (2011). This study focused on trees generated only under a Yule process with a range of λ values. The authors simulated sequences under a simple model of sequence evolution (HKY85), and trees were estimated using BEAST assuming a strict clock and narrow prior or range of prior widths on the root age. Their study assessed the impact of sequence length and nodal calibrations on estimating posterior distributions of λ, and they found that increasing sequence length leads, as expected, to narrower 95% highest posterior density credible intervals of speciation rates. Additionally, broader calibration priors were shown to increase posterior widths of these estimates. It is plausible that forcing estimation of a tree under a particular branching process (such as a Yule process) may impact estimates if the true generating process was different (such as a BD process); this could systematically affect diversification parameter estimates.

Since branch lengths play an important part in estimating diversification parameters, it is also the case that a mismatch of clock models could similarly affect results. Whereas previous work describes a relationship between parameter estimation and misspecification of the model of nucleotide sequence evolution during phylogenetic estimation (Revell, Harmon & Glor, 2005), as well as sequence length and nodal calibrations (Wertheim & Sanderson, 2011), no studies to our knowledge have directly focused on the impact of tree priors and choice of molecular clocks combined (but see Condamine et al., 2015 for comparisons among Yule and BD priors using an empirical dataset). Additionally, a recent study by Duchêne, Hua & Bromham (2017) emphasizes the importance of appropriately accommodating among-lineage molecular rate variation when inferring diversification rates, both of which may be correlated through underlying evolutionary processes. This study simulated datasets under a variety of diversification conditions with a constant background extinction rate and stressed the importance of accurately capturing variable substitution rates as part of reconstructing the phylogeny.

Following a Bayesian statistical philosophy, ideally priors should be selected which reflect a priori knowledge about the data being explored. However, such knowledge may not always be available for each study of interest, especially in non-model systems. It may be possible to use reasonable defaults as selected by an application of choice; however, there is no guarantee that results will always be accurate. One way to tackle this may be to select uninformative (i.e., broad) priors under the assumption that there is enough signal in the data to produce reliable estimates. This can be assessed by performing a parallel analysis sampling only from prior distributions and comparing results to real data or through posterior predictive simulations (Gelman et al., 1995; Huelsenbeck et al., 2001). An alternative approach may be to select the most parameter-rich models with the hope that more complex patterns in the data will be captured and modeled appropriately. However, as the number of parameters increases, issues could arise with overfitting and identifiability, necessitating the use of model selection for comparison of fit. Here, this study is motivated by the observation that phylogenetics, sensu lato, is complex, and the uninitiated may resort to using defaults assigned in tutorials or documentation. At least one other study has explicitly mentioned this (Condamine et al., 2015), referencing selection of a Yule tree prior as suggested in an early BEAST tutorial.

In light of these concerns, we are interested in exploring the choice of tree and molecular clock priors as part of a simulation study conducted using BEAST with choices that researchers may naturally select when interrogating their data. Since we already know that misestimation of the absolute root age of the tree can have dramatic effects on rate estimates, we focus specifically on the effect of priors on relative branch lengths of trees. We quantify the effect of tree prior and clock misspecification on subsequent parameter estimates for diversification models. To accomplish this, we simulate phylogenetic trees and sequence data under a range of combinations of tree priors and molecular clock models. We then re-estimate trees and use these reconstructed trees to calculate maximum likelihood estimates of diversification rate parameters. We compare these estimates to ones from the original trees to evaluate whether or not priors and clock models contribute to error in estimating diversification rates.

Materials and Methods

We take advantage of existing applications to simulate trees under a variety of conditions, simulate nucleotide sequence data on these trees, estimate a tree from the nucleotide data, and estimate diversification rates. The workflow is illustrated in Fig. 1. All scripts are written in the R programming language (R Development Core Team, 2018) and are available on GitHub (https://github.com/bricesarver/prior_simulation_study).

Results

When original trees were simulated under a Yule process, all combinations of tree priors and clocks produced extremely similar estimates to the parameters estimated from trees on which data were simulated (Fig. 2). Distributions overlapped across all combinations of tree priors and molecular clocks. Slight deviations from simulated values are likely attributable to sampling error. The estimates of λ and r were consistently underestimated for the 25-taxa UCLN cases, providing evidence that the number of taxa is important when among-lineage rate heterogeneity is concerned. However, other preliminary trials did not show a consistent pattern of underestimation, suggesting that this pattern results from the 10 trees initially selected for simulation and not a systematic bias. LTT plots of maximum clade credibility trees indicated that the estimated trees generally coincide with the original trees, though the Yule:UCLN case showed greater discordance at nodes deeper in the tree for a small number of replicates (Fig. S1). This is not surprising given the difficulty of estimating nodes deep in the tree, and it also helps explain the discrepancy described above.

When trees were simulated under a BD process, estimates were also generally concordant with the original trees. Medians were nearly identical among many simulation conditions (Fig. 3), though parameters were underestimated in the UCLN cases. This discrepancy was either reduced or did not appear to be present in cases assuming a strict clock. LTT plots revealed that maximum clade credibility trees were, again, approximately equivalent to the original. There were some exceptions, again in the deep nodes of the trees, though these did not drastically affect parameter estimation (Fig. 3). As in the Yule cases, there were no discernable tendencies for parameter estimates to be consistently over or underestimated relative to the simulated trees in preliminary analyses. However, estimates of λ are biased downward, sometimes drastically. For the 25-taxa cases, λ estimates are close to r, even though they ought to be 2r with ε = 0.5. We hypothesize that estimates of λ should approach 2r as the number of taxa increases. To investigate, we performed additional simulations, as described above, but with 50, 75, and 125 taxa. Estimates of λ increase with the number of taxa but are still reduced (Fig. S3).

With low, medium, and high among-lineage substitution rate heterogeneity, assumptions about molecular rates can have substantial impact on parameter estimates (Fig. 4). With low rate heterogeneity, estimates of λ and r are similar to the original trees, but the discordance increases dramatically as the variance in rates among lineages increases. Trees estimated using an UCLN clock appear to suffer the least, especially when estimated under the simulation conditions (BD:UCLN). This effect is most dramatic in the high rate heterogeneity simulations, where the assumption of a tree-wide constant substitution rate can lead to substantially discordant estimates of both λ and r. Further analysis of each of these simulation conditions indicates a deviation from a strict clock, as evidenced by posterior estimates of the coefficient of variation from BEAST on the simulated datasets (95% HPD, low rate hetereogeneity: [0.157–0.1948]; medium rate heterogeneity: [0.2508–1.1833]; high rate heterogeneity: [0.2402–2.6596]).

Discussion

The goal of this study is to determine the impact the choice of tree prior and molecular clock have on the estimation of diversification rates. We focused our efforts on estimating λ, the rate of lineage accumulation, and r, the net diversification rate, under all combinations of two tree priors (Yule and BD) and two flavors of molecular clocks (strict and UCLN). These parameters were selected for investigation because estimating the relative extinction rate (ε) alone is known to be difficult, and estimates of this parameter have larger uncertainty (Nee et al., 1994). Estimating the net diversification rate still provides insight into the effect of extinction across the phylogeny while facilitating a meaningful comparison among simulation conditions. We found that the combination of tree prior and clock did not substantially impact diversification parameter estimates. Across our simulation conditions, parameters from trees estimated under all combinations of tree priors and clocks were concordant with parameter estimates produced from the trees on which nucleotide data were simulated. However, the fact that estimates of λ are biased downward, sometimes drastically, suggests that estimates of λ may be incorrect when trees are estimated assuming no extinction.

The simulations involving low, medium, and high among-lineage substitution rate heterogeneity revealed that it is possible for the choice of clock to have a substantial impact on parameter estimates (Fig. 4). Trees estimated using an UCLN clock appear to suffer the least, especially when estimated under the simulation conditions (BD:UCLN). This effect is most dramatic in the high rate heterogeneity simulations, where the assumption of a tree-wide constant substitution rate can lead to substantially discordant estimates of both λ and r. At the same time, investigators could easily avoid errors associated with using a strict clock by testing for rate heterogeneity in their sequence data.

The assumption of a single rate of evolution across a tree is often violated and can severely impair phylogenetic estimation (Shavit et al., 2007; Penny, 2013). This study assumed rates with a modest amount of heterogeneity, and it appears that a strict clock produces reasonable results in the face of this violation. In other words, a dataset with a small to moderate amount of heterogeneity may have rates that are reasonably captured by a single, global rate. However, it may not be known a priori whether a dataset has disparate rates of evolution among lineages. It would be advisable, then, to assume a clock model that has the potential to model heterogeneity more accurately, and this is partially why the UCLN relaxed clock has seen such widespread use and success in systematic analyses (Drummond et al., 2006). Furthermore, should rates of evolution be extreme among some lineages, it would make sense to attempt to capture any heterogeneity using appropriate priors as opposed to assuming it is absent. Rate homogeneity among lineages, or the absence of a clock altogether, may represent a poor prior given our current understanding of molecular biological processes (Drummond et al., 2006).

There are several caveats to this simulation study. First, our original trees are fully resolved, and nucleotide sequence data are simulated under parameters estimated from a quickly evolving nuclear intron. This implies that there will be a large number of phylogenetically informative sites per individual. Therefore, these trees will be easier to estimate than those that lack signal and/or contain unresolved nodes. Second, there is no extreme rate heterogeneity among lineages. Third, the datasets only contain 25 and 100 taxa, each with only 5,000 bp of nucleotide sequence data, following the protocol of Wertheim & Sanderson (2011). Datasets of this size are considered modest in the current era of high-throughput sequencing, where the generation of hundreds of thousands or millions of base pairs of sequence per sample is possible. More sequence data can lead to more accurate phylogenies, which improves parameter estimates at the expense of computational speed. It is also reasonable to assume that some systems may be best explained through more complex models, such as models that specifically assume multiple, independent diversification rates across a dataset (Alfaro et al., 2009; Rabosky, 2014). Our analyses only assume a single rate of diversification, and this assumption may be violated in larger datasets with greater levels of taxonomic divergence. Furthermore, there are families of models that assume shifts in speciation rates across phylogenies which could be considered (Steel & McKenzie, 2001). Such models can be fit to identify diversification rate heterogeneity and, therefore, estimate diversification rates more accurately under complex scenarios. Finally, by fixing root age, we control for known sources of estimation error that have to do with calibrating molecular evolution when reconstructing time trees. Careful attention to calibrations is essential to obtaining diversification rates in units that are meaningful. We reinforce that it is important to select among models in order to produce accurate, interpretable results for each dataset.

Conclusions

It appears that reasonable parameter estimates can often be achieved regardless of the prior used for phylogenetic tree shape. Among the cases that we simulated, either choice of tree prior appears to capture the underlying branching process; the same holds for molecular clocks with low among-lineage rate heterogeneity. Even in cases with excessive among-lineage rate heterogeneity, it is generally true that existing methods are able to detect and account for that rate variation. Overall, we find that the choice of tree prior and molecular clock has relatively little impact on the estimation of diversification rates.

Supplemental Information