# 12.3: Variation in Diversification Rates through Time

In addition to considering rate variation across clades, we might also wonder whether birth and/or death rates have changed through time. For example, perhaps we think our clade is an adaptive radiation that experienced rapid diversification upon arrival to an island archipelago and slowed as this new adaptive zone got filled (Schluter 2000). This hypothesis is an example of density-dependent diversification, where diversification rate depends on the number of lineages that are present (Rabosky 2013). Alternatively, perhaps our clade has been experiencing extinction rates that have changed through time, perhaps peaking during some time period of unfavorable climactic conditions (Benton 2009). This is another hypothesis that predicts variation in diversification (speciation and extinction) rates through time.

We can fit time-dependent diversification models using likelihood equations that allow arbitrary variation in speciation and/or extinction rates, either as a function of time or depending on the number of other lineages in the clade. To figure out the likelihood we can first make a simplifying assumption: though diversification rates might change, they are constant across all lineages at any particular time point. In particular, this means that speciation (and/or extinction) rates slow down (or speed up) in exactly the same way across all lineages in an evolving clade. This is again the Equal-Rates Markov (ERM) model for tree growth described in the previous chapter.

Our assumption about equal rates across lineages at any time means that we can consider time-slices through the tree rather than individual branches, i.e. we can get all the information that we need to fit these models from lineage through time plots.

The most general way to fit time-varying birth-death models to phylogenetic trees is described in Morlon et al. (2011). Consider the case where both speciation and extinction rates vary as a function of time, λ(t) and μ(t). Morlon et al. (2011) derive the likelihood for such a model as:

$$L(t_1, t_2, \dots, t_{n-1}) = (n+1)! \dfrac{f^n \sum_{i=1}^{n-1} \lambda(t_i) \Psi(s_{i,1}, t_i) \Psi(s_{i,2}, t_i)}{\lambda [1-E(t_1)]^2} \label{12.1}$$

where:

$$E(t) = 1 - \frac{e^{\int_0^t [\lambda(u) - \mu(u)]du}}{\dfrac{1}{f} + \displaystyle \int_0^t e^{\int_0^s [\lambda(u) - \mu(u)]du} \lambda(s) ds}\label{12.2}$$

and:

$$\Psi(s, t) = e^{\int_s^t [\lambda(u) - \mu(u)]du} \left[ 1 + \frac{\displaystyle \int_s^t e^{\int_0^\tau [\lambda(\sigma) - \mu(\sigma)]d\sigma}\lambda(\tau)d\tau}{\frac{1}{f} + \displaystyle \int_0^s e^{\int_0^\tau \left[\lambda(\sigma) - \mu(\sigma) \right]d\sigma} \lambda(\tau) d\tau} \right]^{-2} \label{12.3}$$

Following Chapter 11, n is the number of tips in the tree, and divergence times t1, t2, …, tn − 1 are defined as measured from the present (e.g. decreasing towards the present day). λ(t) and μ(t) are speciation and extinction rates expressed as an arbitrary function of time, f is the sampling fraction (under a uniform sampling model). For a node starting at time ti, si, 1 and si, 2 are the times when the two daughter lineages encounter a speciation event in the reconstructed tree1. E(t), as before, is the probability that a lineage alive at time t leaves no descendants in the sample. Finally, $$Ψ(s, t)$$ is the probability that a lineage alive at time s leaves exactly one descent at time t < s in the reconstructed tree. These equations look complex, and they are - but basically involve integrating the speciation and extinction functions (and their difference) along the branches of the phylogenetic tree.

Note that my equations here differ from the originals in Morlon et al. (2011) in two ways. First, Morlon et al. (2011) assumed that we have information about the stem lineage and, thus, used an index on ti that goes up to n instead of n − 1 and a different denominator conditioning on survival of the descendants of the single stem lineage (Morlon et al. 2011). Second, I also multiply by the total number of topological arrangements of n taxa, (n + 1)!.

If one substitutes constants for speciation and extinction (λ(t)=λc, μ(t)=μc) in Equation \ref{12.1}, then one obtains equation 11.24; if one additionally considers the case of complete sampling and substitutes f = 1 then we obtain equation 11.18. This provides a single unified framework for time-varying phylogenetic trees with uniform incomplete sampling (see also Höhna 2014 for independent but equivalent derivations that also extend to the case of representative sampling).

Equation \ref{12.1} requires that we define speciation rate as a function of time. Two types of time-varying models are currently common in the comparative literature: linear and expoential. If speciation rates change linearly through time (see Rabosky and Lovette 2008 for an early version of this model):

$λ(t)=λ_0 + α_λt \label{12.4}$

Where λ0 is the initial speciation rate at the present and alpha is the slope of speciation rate change as we go back through time. αλ must be chosen so that speciation rates do not become negative as we move back through the tree: αλ > −λ0/t1. Note that the interpretation of αλ is a bit strange since we measure time backwards: a positive αλ, for example, would mean that speciation rates have declined from the past to the present. Other time-dependent models published earlier (e.g. Rabosky and Lovette 2008, which considered a linearly declining pure-birth model) do not have this property.

We could also consider a linear change in extinction through time:

$μ(t)=μ_0 + α_μt \label{12.5}$

Again, αμ is the change in extinction rate through time, and must be interpreted in the same "backwards" way as αλ. Again, we must restrict our parameter to avoid a negative rate: αμ > −μ0/t1

One can then substitute either of these formulas into Equation \ref{12.1} to calculate the likelihood of a model where speciation rate declines through time. Many implementations of this approach use numerical approximations rather than analytic solutions (see, e.g., Morlon et al. 2011; Etienne et al. 2012).

Another common model has speciation and/or extinction rates changing exponentially through time:

$λ(t)=λ_0e^{β_λt} \label{12.6}$

and/or

$μ(t)=μ_0e^{β_μt} \label{12.7}$

We can again calculate likelihoods for this model numerically (Morlon et al. 2011, Etienne et al. (2012)).

As an example, we can test models of time-varying diversification rates across part of the amphibian tree of life from Jetz and Pyron (2018). I will focus on one section of the salamanders, the lungless salamanders (Plethodontidae, comprised of the clade that spans Bolitoglossinae, Spelerpinae, Hemidactylinae, and Plethodontinae). This interesting clade was already identified above as including both an increase in diversification rates (at the base of the clade) and a decrease (on the branch leading to Hemidactylinae; Figure 12.4). The tree I am using may be missing a few species; this section of the tree includes 440 species in Jetz and Pyron (2018) but 471 speices are listed on Amphibiaweb as of May 2018. Thus, I will assume random sampling with f = 440/471 = 0.934. Comparing the fit of a set of models, we obtain the following results:

Model Number of params. Param. estimates lnL AIC
CRPB 1 λ = 0.05111267 497.8 -993.6
CRBD 2 λ = 0.05111267 497.8 -991.6
SP-L μ = 0
3 λ0 = 0.035 513.0 -1019.9
αλ = 0.0011
μ = 0
SP-E 3 λ0 = 0.040 510.7 -1015.4
βλ = 0.016
μ = 0
EX-L 3 λ = 0.053 497.8 -989.6
μ0 = 0
αμ = 0.000036
EX-E 3 λ = 0.069 510.6 -1015.3
μ0 = 61.9
βμ = −111.0

Models are abbreviated as: CRPB = Constant rate pure birth; CRBD = Constant rate birth-death; SP-L = Linear change in speciation; SP-E = Exponential change in speciation; EX-L = Linear change in extinction; EX-E = Exponential change in extinction.

The model with the lowest AIC score has a linear decline in speciation rates, and moderate support compared to all other models. From this, we support the inference that diversification rates among these salamanders has slowed through time. Of course, there are other models I could have tried, such as models where both speciation and extinction rates are changing through time, or models where there are many more extant species of salamanders than currently recognized. The conclusion we make is only as good as the set of models being considered, and one should carefully consider any plausible models that are not in the candidate set.