# Ranked Tree Shapes, Nonrandom Extinctions, and the Loss of Phylogenetic Diversity

Ranked Tree Shapes, Nonrandom Extinctions, and the Loss of Phylogenetic Diversity Abstract Phylogenetic diversity (PD) is a measure of the evolutionary legacy of a group of species, which can be used to define conservation priorities. It has been shown that an important loss of species diversity can sometimes lead to a much less important loss of PD, depending on the topology of the species tree and on the distribution of its branch lengths. However, the rate of decrease of PD strongly depends on the relative depths of the nodes in the tree and on the order in which species become extinct. We introduce a new, sampling-consistent, three-parameter model generating random trees with covarying topology, clades relative depths, and clades relative extinction risks. This model can be seen as an extension to Aldous’ one parameter splitting model ($$\beta$$, which controls for tree balance) with two additional parameters: a new parameter $$\alpha$$ quantifying the relation between age and richness of subclades, and a parameter $$\eta$$ quantifying the relation between relative abundance and richness of subclades, taken herein as a proxy for overall extinction risk. We show on simulated phylogenies that loss of PD depends on the combined effect of all three parameters, $$\beta$$, $$\alpha,$$ and $$\eta$$. In particular, PD may decrease as fast as species diversity when high extinction risks are clustered within small, old clades, corresponding to a parameter range that we term the “danger zone” ($$\beta<-1$$ or $$\alpha<0$$; $$\eta>1$$). Besides, when high extinction risks are clustered within large clades, the loss of PD can be higher in trees that are more balanced ($$\beta>0$$), in contrast to the predictions of earlier studies based on simpler models. We propose a Monte-Carlo algorithm, tested on simulated data, to infer all three parameters. Applying it to a real data set comprising 120 bird clades (class Aves) with known range sizes, we show that parameter estimates precisely fall close to the danger zone: the combination of their ranking tree shape and nonrandom extinctions risks makes them prone to a sudden collapse of PD. Beta-splitting model, biodiversity, broken stick, field of bullets model, macroevolution, phylogenetic tree, rarefaction, sampling distribution, self-similar fragmentation As it becomes increasingly clear that human activities are causing a major extinction crisis (Leakey and Lewin, 1995; Glavin, 2007; Wake and Vredenburg, 2008; Barnosky et al., 2011), several theoretical studies have aimed at characterizing how the evolutionary legacy of parts of the tree-of-life, and hence also the genetic diversity able to drive future evolution, will decrease in the face of forthcoming extinctions. This evolutionary component of biodiversity can be measured by the phylogenetic diversity (PD), defined as the sum of the branch lengths of the phylogeny spanned by a given set of taxa (Faith, 1992). This metric is increasingly being used to measure biodiversity and to identify conservation strategies (Veron et al., 2015). Nee and May (1997) were the first to formally investigate the expected loss of PD in the face of species extinctions, by simulating species trees using the Kingman coalescent. They found that 80% of the PD can be conserved even when 95% of species are lost. Further studies showed that the loss of PD is in fact much higher when trees are generated through other models of species diversification, such as the Yule or the birth–death models (Morlon et al., 2011b; Mooers et al., 2012; Lambert and Steel, 2013). These models indeed generate longer pendant edges (i.e., branches that lead to the tips), hence lower phylogenetic redundancy, than in the standard Kingman coalescent (used by Nee and May, 1997). However, Nee and May (1997) also showed that PD is very sensitive to the shape of the species tree (also called its “topology”), with extremely unbalanced trees (“caterpilar trees”) losing much more PD than balanced trees (“bush trees”), due to a lack of phylogenetic redundancy (i.e., the presence of recently diverged sister species). Overall, these results highlighted the sensitivity of the loss of PD in response to species extinctions to both edge lengths and tree shape. In this line, we also expect the correlation between the species richness of clades and their relative ages to have a significant impact on the loss of PD (“clade” standing here for any subtree within the full phylogeny). Here the age of a clade, also called “stem age,” denotes the depth (measured from the present) of its root node (i.e., the node where this clade is tied to the rest of the tree). Under random extinction, since smaller clades are more likely to become extinct first, the consequence of their total extinction on PD will depend on the lengths of pendant edges in these clades compared to those in larger clades. Under models with diversification rates that are constant through time and homogeneous across lineages, the time for speciation hypothesis states that the size of clades is correlated to their age (Magallon and Sanderson, 2001), yet several empirical studies found no correlation between the two (Ricklefs 2007; Rabosky et al. 2012; Sánchez-Reyes et al. 2016; but see McPeek and Brown 2007). The effect of such correlation on the loss of PD has not yet been explored, but should be particularly important in unbalanced phylogenetic trees (exhibiting large variation in the species richness of clades), which dominate empirical data (e.g., Guyer and Slowinski, 1991; Heard, 1992; Guyer and Slowinski, 1993; Slowinski and Guyer, 1993; Mooers, 1995; Purvis, 1996; Mooers and Heard, 1997; Blum and François, 2006). Besides, the loss of PD was shown to be influenced by the distribution of extinction risks within species trees. Several studies showed that accounting for realistic scenarios of species extinctions (considering that species with higher extinction risk as per the International Union for Conservation of Nature (IUCN) Red List status are more likely to go extinct first) predicts proportionately higher losses in PD than scenarios with random extinction risks (e.g., and review, Purvis et al., 2000a; von Euler, 2001; Purvis, 2008; Veron et al., 2015). Extinctions may for example be clustered within certain clades (Bennett and Owens, 1997; McKinney, 1997; Russell et al., 1998; Purvis et al., 2000a; Baillie et al., 2004; Bielby et al., 2006; Fritz and Purvis, 2010), correlated to the age of clades (von Euler, 2001; Johnson et al., 2002; Redding and Mooers, 2006), to the species richness of clades (Russell et al., 1998; Hughes, 1999; Purvis et al., 2000a; Schwartz and Simberloff, 2001; von Euler, 2001; Johnson et al., 2002; Lozano and Schwartz, 2005; Vamosi and Wilson, 2008, assuming in some studies a correlation between rarity and extinction risks), or to speciation rates (Heard and Mooers, 2000). In contrast, most of the theoretical analyses of predictions based on model trees (Nee and May, 1997; Mooers et al., 2012; Lambert and Steel, 2013) have all been based so far on the field of bullets model, which considers equal extinction probabilities across species (Raup et al., 1973; Van Valen, 1976; Nee and May, 1997; Vazquez and Gittleman, 1998; but see Heard and Mooers 2000). One can assume extinction events are independent but not identically distributed across species, as considered in the generalized field of bullets model (Faller et al., 2008). In an exchangeable phylogenetic model in which extinction probabilities are themselves random and independent with the same distribution, this would not affect the overall loss of PD (as both models are stochastically equivalent, Lambert and Steel, 2013). However, as stated by Faller et al. (2008), it is essential to explore models that weaken the strong assumption in the (generalized) field of bullets models that extinction events are randomly and independently distributed among the tips of phylogenetic trees. Here, we hence investigate how the loss of PD is influenced by the two abovementioned factors: (i) the ranked shape of the species tree, characterized by the relation between the richness of clades and their age or depth in the tree and (ii) nonrandom extinctions, characterized by the relation between the richness of clades and the extinction risks within them. Here, “ranked shape” refers to the shape of the tree combined with the additional knowledge of relative depths—the order in which nodes appear in the tree, but to the exclusion of the actual divergence times (e.g., Lambert et al., 2017). We introduce a three-parameter model generating random ranked tree shapes endowed with random numbers summing to one at the tips, interpreted as relative abundances (or geographic ranges) of contemporary species. This model can be seen as an extension to Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001) with two additional parameters: a parameter $$\alpha$$ quantifying the relation between the richness of a clade and its relative age (i.e., the rank of appearance of its root node) termed “age-richness index” hereafter, and another parameter $$\eta$$ quantifying the relation between the richness of a clade and its relative abundance or frequency (i.e., the sum of abundances of the species it encompasses divided by the sum of abundances of all extant species in the phylogeny), termed “abundance-richness index” hereafter. When $$\beta=0$$ and $$\alpha=1$$, the ranked shape of the tree is the same as the ranked shape of a standard coalescent tree or of a Yule tree stopped at a fixed time (see Proposition 1 in Supplementary Appendix S1 available on Dryad). We further assume that extinctions of contemporary species occur sequentially in the order of their abundances, starting with the least abundant species, which roughly reduces to the field of bullets model when $$\eta=1$$ (see Proposition 2 in Supplementary Appendix S1 available on Dryad). The parameters of the model are not supposed to map onto biological processes. Our aim is to produce and describe a broad range of ranked tree shapes and extinction risk distributions. The model nevertheless reflects relevant biological patterns. The imbalance and node order of phylogenies may be affected by diversification rates varying with time (Rabosky and Lovette, 2008b; Moen and Morlon, 2014), with species age (Doran et al., 2006; Hagen et al., 2015; Alexander et al., 2015), or among lineages (Cardillo et al., 2005; Maddison et al., 2007; Alfaro et al., 2009; Morlon et al., 2011a). Extinction risk may cluster within the phylogeny if it is correlated with some species characteristic, such as body size (Gaston and Blackburn, 1995; Johnson et al., 2002) or habitat use (Johnson et al., 2002). Finally, tree imbalance, node order and extinction risk clustering are likely to interact if differences in diversification rates and extinction risks are driven by the same trait, or are the result of a common process (e.g., geographic speciation, Pigot et al. 2010). We explore the rate of decrease of PD as species sequentially become extinct, based on simulated data under variation in all three parameters over a significant range of their possible values. Interestingly, the joint variation of the parameter $$\eta$$ with the ranked shape of species trees (set by parameters $$\beta$$ and $$\alpha$$) affects the clustering of extinction risks and the relationship between extinction risks and clade age (determined by the similarity or dissimilarity of the direction of deviations of $$\alpha$$ and $$\eta$$ from 1). Therefore, considering simultaneous variation in $$\beta$$, $$\alpha,$$ and $$\eta$$ allows us to explore the effects on the loss of PD of the different patterns of nonrandom extinctions observed in empirical data. We therefore provide general predictions on the sensitivity of the evolutionary legacy of clades to extinction, as a function of three simple statistics summarizing tree balance, ranked tree shape and the distribution of extinction risks across clades. Following this exploration, we then propose a Monte-Carlo inference algorithm enabling maximum likelihood estimation of the parameters $$\beta$$, $$\alpha,$$ and $$\eta$$ from real data sets. When tested against simulated data, this algorithm performs reasonably well over a wide range of parameter values for phylogenies with 50 tips or more. The estimates of parameters (beta, alpha, eta) on a real data set of bird family phylogenies and their range size distributions finally reveal empirical patterns clustered within a given parameter zone which make these clades particularly prone to strong loss of PD. Materials and Methods Modeling Ranked Tree Shapes The first version of the model we present allows one to generate random ranked tree shapes, that is tree shapes endowed with the additional knowledge of node ranks. Usually, one can generate random ranked tree shapes by time-continuous branching processes stopped at some fixed or random time, where particles are endowed with a heritable trait influencing birth and death rates. In these models, it is generally not possible to characterize the distribution of the tree shape (for an exception, see Sainudiin and Véber, 2016) or to relate it to known distributions whenever it does not have the shape of the Yule tree (i.e., the tree generated by a pure-birth process). Also, since the same trait is usually responsible for both the tree shape and the order of nodes, it is impossible to disentangle the roles of either of these characteristics on the behavior of the tree in the face of current extinctions. Last, these models do not fulfill a very important property called sampling consistency (usually considered in combination with exchangeability, i.e., ecological equivalence between species). This property ensures that one can equivalently draw a random tree with $$n$$ tips from the distribution or draw a tree with $$n+1$$ tips and then remove one tip at random. The model we propose here has two parameters: $$\beta \in (-2,+\infty)$$ (tree balance index) determines the balance of the tree, similarly as in Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001), and $$\alpha \in (-\infty,+\infty)$$ (age-richness index) sets the relation between the species richness of a clade and its relative age (Fig. 2). Figure 1. View largeDownload slide Illustration of the model generating ranked tree shapes. Construction of the ranked shape of a tree containing $$N=5$$ species. (1) Five random marks $$(U_i)_{i \in \lbrace 1,...,5\rbrace}$$ are drawn uniformly in the interval $$[0,1]$$ (the marks on the bottom line). (2) At each time step (time flowing downwards), we randomly select one interval $$X$$, with each interval $$X_j$$ having a weight $$|X_j|^{\alpha}$$ (in black). Then, we draw a random variable $$R$$ in a beta distribution with parameters ($$\beta+1$$, $$\beta+1$$), and split the selected interval $$X$$ into two subintervals, $$X_{\rm{left}}$$ of size $$R|X|$$ and $$X_{\rm{right}}$$ of size $$(1-R)|X|$$ (dark red mark). (3) Repeating this process over time until all intervals $$X_j$$ contain only one mark leads a tree with a ranked shape. Dotted branches correspond to unsampled subtrees (i.e., there is no mark in the corresponding interval). Figure 1. View largeDownload slide Illustration of the model generating ranked tree shapes. Construction of the ranked shape of a tree containing $$N=5$$ species. (1) Five random marks $$(U_i)_{i \in \lbrace 1,...,5\rbrace}$$ are drawn uniformly in the interval $$[0,1]$$ (the marks on the bottom line). (2) At each time step (time flowing downwards), we randomly select one interval $$X$$, with each interval $$X_j$$ having a weight $$|X_j|^{\alpha}$$ (in black). Then, we draw a random variable $$R$$ in a beta distribution with parameters ($$\beta+1$$, $$\beta+1$$), and split the selected interval $$X$$ into two subintervals, $$X_{\rm{left}}$$ of size $$R|X|$$ and $$X_{\rm{right}}$$ of size $$(1-R)|X|$$ (dark red mark). (3) Repeating this process over time until all intervals $$X_j$$ contain only one mark leads a tree with a ranked shape. Dotted branches correspond to unsampled subtrees (i.e., there is no mark in the corresponding interval). Figure 2. View largeDownload slide Phylogenetic trees simulated for different values of $$\beta$$ (tree balance) and $$\alpha$$ (age-richness index). Node depths are set as in a Yule pure-birth process. Parameter values: $$\beta=-1.5$$ (bottom) or 10 (top), $$\alpha=-10$$ (left) or 10 (right), number of species $$N=30$$, approximation parameter (see Supplementary Appendix S1 available on Dryad) $$\varepsilon=0.001$$. Figure 2. View largeDownload slide Phylogenetic trees simulated for different values of $$\beta$$ (tree balance) and $$\alpha$$ (age-richness index). Node depths are set as in a Yule pure-birth process. Parameter values: $$\beta=-1.5$$ (bottom) or 10 (top), $$\alpha=-10$$ (left) or 10 (right), number of species $$N=30$$, approximation parameter (see Supplementary Appendix S1 available on Dryad) $$\varepsilon=0.001$$. The construction of a tree according to this model is done by following the steps indicated hereunder (illustrated on Fig. 1). We start with $$n$$ uniform, independent random variables $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$ in the interval $$[0,1]$$. For each $$i$$, the mark $$U_i$$ is associated to the tip species labeled $$i$$ in the phylogeny. The procedure consists in sequentially partitioning $$[0,1]$$ into a finite subdivision thanks to random variables independent of the marks $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$, until all marks are in distinct components of the partition. At each step, the new point added to the subdivision corresponds to a split event in the tree. In the beginning, there is only one component in the partition (the interval $$[0,1]$$ itself). Each interval $$X$$ of the partition containing at least two marks among the $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$ is given a weight equal to $$|X|^\alpha$$, where $$|X|$$ denotes the width of the interval $$X$$. One of these intervals is selected with a probability proportional to its weight. Draw a random variable $$R$$ in a Beta distribution with parameters $$(\beta+1, \beta+1)$$. The selected interval $$X$$ of width $$|X|$$ is then split into two disjoint subintervals, $$X_{\rm{left}}$$ and $$X_{\rm{right}}$$, with widths $$|X_{\rm{left}}|=R|X|$$ and $$|X_{\rm{right}}|=(1-R)|X|$$. Each subinterval contains a distinct subset of the marks. The marks in the subinterval $$X_{\rm{left}}$$ determine the tips in the left subtree of the phylogeny, and the marks in the subinterval $$X_{\rm{right}}$$ determine the tips in the right subtree. This step is performed even if one subinterval contains no mark among the $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$, which corresponds to a subtree with no sampled species (i.e., species which are seen in the phylogeny). The order in which the splitting subintervals are selected sets the order of branching events (i.e., nodes) in the tree. If no interval contains more than one mark, the process is stopped. Otherwise, go to Step 1. We can relate the tree shape in this model to well-known distributions. Because $$\alpha$$ has no impact on the way we refine the subdivision, the tree shape generated with our model coincides exactly with the tree shape obtained in Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001) with the same parameter $$\beta$$. For small values of $$\beta$$, the intervals are often split close to one of their extremities, and the resulting tree is unbalanced, converging to the perfectly unbalanced “caterpilar” tree as $$\beta\to -2$$. On the contrary, for large values of $$\beta$$, the intervals are often split close to their middle, and the resulting tree is balanced. We stress that unlike most models, $$\alpha$$ can be tuned independently of $$\beta$$, allowing node ranks to vary while keeping the same tree shape. For small values of $$\alpha$$ (in particular $$\alpha<0$$), the smallest subintervals have a higher probability of being selected, so smaller clades tend to be older. On the contrary, for large values of $$\alpha$$, the largest subintervals have a higher probability of being selected, so smaller clades tend to be younger. We notice that as $$\beta$$ gets close to $$-2$$ the effect of $$\alpha$$ vanishes, since most of the time there is only one subinterval containing more than one mark, so only one subinterval to split. In maximally unbalanced tree shape ($$\beta=-2$$), there is only one ranked tree shape and the order of nodes is fixed, so $$\alpha$$ plays no role. As is well-known, the tree obtained with $$\beta=0$$ has the same shape has the tree generated with the Yule process (Yule, 1925) or the Kingman coalescent (Kingman, 1982) after ignoring node ranks (Nee, 2006; Lambert and Stadler, 2013). When $$\alpha=1$$ in addition to $$\beta=0$$, we show in Supplementary Appendix S1 (Proposition 1) available on Dryad at http://dx.doi.org/10.5061/dryad.mv980, that our model generates the same tree shape with node ranks as Yule trees, which is actually known to be the same as the ranked tree shape of the Kingman coalescent tree. The version of the model we present here only allows simulation of trees with $$\beta>-1$$, as the beta distribution is only defined for positive parameter values. Actually, our model coincides with the ranked tree in a self-similar, binary fragmentation with self-similarity index $$\alpha$$ (which is the very age-richness index $$\alpha$$) and with fragmentation measure $$\int_0^1 \delta_{(x,1-x,0,0,\ldots)}\,x^{\beta+1}(1-x)^{\beta+1}dx$$ (as defined in Bertoin, 2002, 2006), which makes sense as soon as $$\beta>-2$$. In Supplementary Appendix S1 (Proposition 3) available on Dryad, we present an algorithm based on fragmentation processes equivalent to that presented above (using one additional approximation parameter $$\varepsilon$$, the maximal frequency of unsampled clades with insignificant richness, consistently set to 0.001). Albeit less intuitive, this method allows us to simulate trees for all $$\beta>-2$$. Last, it is important to notice that our model is both exchangeable and sampling consistent. It is exchangeable because labels can be swapped without changing the distribution of the tree, since marks all have the same distribution. It is sampling consistent because removing tip labeled $$n+1$$ (or any other tip, by exchangeability) amounts to removing mark $$U_{n+1}$$, which does not modify the ranked tree shape obtained from marks $$(U_i)_{i \in \lbrace1,...,n\rbrace}$$. Incorporating Nonrandom Extinctions In order to parametrize the relation between the richness of a clade and its relative abundance, we now add to the model a new parameter $$\eta\ge 0$$ called “abundance-richness index.” Each time an interval $$X$$ is split into two subintervals, $$X_{\rm{left}}$$ and $$X_{\rm{right}}$$ with widths $$|X_{\rm{left}}|=R|X|$$ and $$|X_{\rm{right}}|=(1-R)|X|$$, each of the two subtrees is granted a part of the relative abundance $$A_X$$ of the parental clade equal to   \begin{align*} A_{X_{\rm{left}}} & = \frac{|X_{\rm{left}}|^\eta}{ |X_{\rm{left}}|^\eta+|X_{\rm{right}}|^\eta} \,A_X = \frac{R^\eta}{R^\eta+(1-R)^\eta}\, A_X \\ A_{X_{\rm{right}}} & = \frac{|X_{\rm{right}}|^\eta}{ |X_{\rm{left}}|^\eta+|X_{\rm{right}}|^\eta}\, A_X = \frac{(1-R)^\eta}{R^\eta+(1-R)^\eta} \, A_X. \end{align*} This way of allocating frequencies to taxa is reminiscent of the “broken stick model” (MacArthur, 1957; MacArthur and Wilson, 1967; Colwell and Lees, 2000), where the unit interval is broken into subintervals each representing the frequency or resource share of each species or clade in the community, even though in our case the allocation of relative abundances has no mechanistic interpretation. This is usually done by throwing uniform points independently in the interval or by throwing the points sequentially, always to the right of the last one, leading to the Poisson–Dirichlet distribution appearing in mathematical population genetics (Feng, 2010; Ewens, 2012) as well as in the neutral theory of biodiversity (Hubbell, 2001). The model remains sampling-consistent insofar as each $$A_X$$ is interpreted as the relative abundance of an entire clade, that is the sum of relative abundances of all species, sampled or not (seen or not in the phylogeny), belonging to this clade. Sampling consistency now means that generating a ranked tree shape with relative abundances on $$n$$ tips is equivalent to the following process: generate a ranked tree shape with relative abundances on $$n+1$$ tips, remove one tip at random and sum the abundance of the removed tip to that of its sister clade (i.e., the clade descending from the interior node connected to the removed tip by a pendant edge). In the extinction numerical experiment, we determine the order of species extinctions deterministically based on their abundance: rarest species become extinct first, whereas most abundant species become extinct last. The distribution of relative abundances across species is captured in the abundance-richness index $$\eta$$ (Fig. 3 and in the Supplementary Appendix S1 available on Dryad): when $$\eta=1$$, all tip species have the same relative abundance in expectation ($$\frac{1}{n}$$), such that the correlation between the relative abundance of a clade and its richness would also be close to $$1.0$$. When we simulate extinction with $$\eta=1$$, all species have an approximately equal chance of being removed, and we approach a field of bullets model (see Proposition 2 in Supplementary Appendix S1 available on Dryad; in the case $$\beta=-1$$, the equivalence is exact). For $$\eta>1$$, species in larger clades have larger abundances on average (and thus lower extinction risks), and for $$\eta<1$$, species in larger clades have smaller abundances (and thus higher extinction risks) on average. This modeling approach allows us to tune the sign and strength of the relation between the richness of a clade and the extinction risk of its tip species. Figure 3. View largeDownload slide Distribution of species frequencies across the tips of phylogenetic trees for different values of abundance-richness index $$\eta$$. Dot sizes sort species according to their frequency (larger dots for more abundant species). Parameter values: $$\eta=0.2$$, 1, or 3 (from left to right), $$\beta=0$$, $$\alpha=0$$, number of species $$N=30$$, approximation parameter $$\varepsilon=0.001$$. Results with $$\beta=-1.9$$ are shown in the Supplementary Appendix S1 available on Dryad. Figure 3. View largeDownload slide Distribution of species frequencies across the tips of phylogenetic trees for different values of abundance-richness index $$\eta$$. Dot sizes sort species according to their frequency (larger dots for more abundant species). Parameter values: $$\eta=0.2$$, 1, or 3 (from left to right), $$\beta=0$$, $$\alpha=0$$, number of species $$N=30$$, approximation parameter $$\varepsilon=0.001$$. Results with $$\beta=-1.9$$ are shown in the Supplementary Appendix S1 available on Dryad. Testing the Effect of $$\beta$$, $$\alpha$$, and $$\eta$$ on PD Loss The effect of all three model parameters on the relationship between species loss and PD loss is studied in a systematic way by simulation. We considered values of $$\beta$$ in $$(-2,10]$$, values of $$\alpha$$ in $$[-3,3]$$ and $$\eta$$ in $$[0.1,3]$$. Because our model specifies how interior nodes are ranked in time but not their actual timing, we use a pure-birth process to generate node depths, adding the latter on top of ranked tree shapes. The use of another model for generating node depths leads to qualitatively similar results, albeit quantitatively different (as an illustration, we show results with edge lengths set as in the Kingman coalescent in Supplementary Appendices S4 and S6 available on Dryad). For each set of parameter values, we generated 100 trees with 100 tips ($$N=100$$). We sequentially removed extinct species from these trees (in the order of increasing species abundances, as explained earlier), and computed the remaining PD (sum of all branch lengths; Faith, 1992) for increasing fractions of extinct species. Parameter Inference We inferred the parameters $$\beta$$, $$\alpha,$$ and/or $$\eta$$ from simulated or empirical data sets by maximum likelihood. As is already well-known (Aldous, 1996; Blum and François, 2006; Lambert et al., 2017), the likelihood of a labeled tree shape under Aldous’ $$\beta$$-splitting model is explicit. Since the likelihood of the tree shape under our model is the same as in Aldous’ model (and in particular independent of $$\alpha$$ and $$\eta$$), we can use it to estimate $$\beta$$. In contrast, computing the likelihood of the ranked tree shape requires to follow through time the lengths of all intervals of the partition containing marks, which may decrease without separating marks (unsampled species). Given that the likelihood of the ranked tree (with or without tip abundances) with the additional knowledge of interval lengths is explicit, we use a Monte-Carlo data augmentation procedure, in which the augmentation data are the numbers and sizes of unsampled splits on each branch (which allow us to reconstruct the interval lengths through time). The likelihood of the ranked tree with tip abundances is then computed by averaging over augmentations and is optimized over possible values of $$(\alpha, \eta)$$. We first tested our ability to infer the model parameters on simulated trees. To do so we simulated trees with 20, 50 and 100 tips for all possible combinations of $$\alpha$$ in $$\lbrace -1 , 0, 1, 2 \rbrace$$, $$\beta$$ in $$\lbrace -1 , 0, 1 \rbrace$$ and $$\eta$$ in $$\lbrace 0.2 , 0.5, 1, 1.5 , 2 \rbrace$$. For each tree size and parameter combination, we simulated 20 trees with tip abundances, for a total number of 3600 trees. We then inferred the model parameters on these trees and compared them to the values used in the simulations. The inference of the parameter $$\beta$$ was straightforward, being computed as the maximum likelihood estimate on the interval $$\left[ -2,10\right]$$ with the function maxlik.betasplit from the R-package apTreeshape (Bortolussi et al., 2006). The parameters $$\alpha$$ and $$\eta$$ were estimated with the method introduced hereabove, with values respectively constrained in the intervals $$\left[-4,4 \right]$$ and $$\left[0.1,10 \right]$$. The value of $$\varepsilon$$ (approximation parameter, see Supplementary Appendix S1 available on Dryad) was here again fixed to $$0.001$$. After validating the estimation procedure, we applied it to real bird family trees. We used the Maximum Clade Credibility (MCC) tree from Jetz et al. (2012), and pruned it to keep family level phylogenies. We kept only the phylogenies that included at least 50 species, and used range sizes from Map of Life (https://mol.org/) as proxies for relative abundances. The value of $$\varepsilon$$ and the constraints on parameter ranges were here the same as in the test on simulated phylogenies. The model was coded—and the analyses of phylogenetic trees were performed—using R (R Development Core Team, 2012) and the R packages cubature (Johnson and Narasimhan, 2013), ape (Paradis et al., 2004), sads (Prado et al., 2015), apTreeshape (Bortolussi et al., 2006), and picante (Kembel et al., 2014). The code is available in the R-package apTreeshape (Bortolussi et al., 2006). The list of available functions is given in Supplementary Appendix S10 available on Dryad. Results Influence of Ranked Tree Shape on PD Loss Here we only address the influence of $$\alpha$$ on PD loss, assuming a field of bullets model for species extinctions ($$\eta=1$$). The expected PD loss is then a convex function of the fraction $$p$$ of extinct species (as proved mathematically for any binary tree under the field of bullets model, see Eq (34) in Lambert and Steel, 2013), always lying below $$p$$ (Fig. 4A,C,E,G). Figure 4. View largeDownload slide Influence of the ranked tree shape (tree balance $$\beta$$ and age-richness index $$\alpha$$) on PD loss, for increasing fractions of species extinctions. Tree balance $$\beta$$ changes from $$10$$ (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”). Results are shown either as a function of the extinction fraction $$p$$ (left column; for different $$\alpha$$ values) or as a function of $$\alpha$$ (right column; for different extinction fractions $$p$$). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right in a, c, e, g; from light to dark color in b, d, f, h). The dotted lines in a, c, e, g show the bisector. Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 4. View largeDownload slide Influence of the ranked tree shape (tree balance $$\beta$$ and age-richness index $$\alpha$$) on PD loss, for increasing fractions of species extinctions. Tree balance $$\beta$$ changes from $$10$$ (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”). Results are shown either as a function of the extinction fraction $$p$$ (left column; for different $$\alpha$$ values) or as a function of $$\alpha$$ (right column; for different extinction fractions $$p$$). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right in a, c, e, g; from light to dark color in b, d, f, h). The dotted lines in a, c, e, g show the bisector. Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Consistently with previous studies (Nee and May, 1997; von Euler, 2001) we find that when the relation between age and richness of clades is similar as that in Yule trees ($$\alpha=1$$), very unbalanced trees (caterpilar-like trees) lose more PD in the face of species extinctions than Yule or more balanced trees (Fig. 4G,H vs. A–D, with $$\alpha=1$$). The effect is nonlinear in $$\beta$$: the tree shape has little influence on the loss of PD when $$\beta\ge -1$$, but increases sharply as $$\beta$$ decreases from $$-1$$ to $$-1.9$$ (results as a function of $$\beta$$ in Supplementary Appendix S2 available on Dryad). Unbalanced tree shapes are associated with the presence of long edges leading to evolutionary distinct species (Fig. 2). These edges constitute an important fraction of the PD in unbalanced species trees, so that their extinction generates a significant drop in PD. As $$\beta$$ gets closer to $$-2$$ (case of the “caterpilar tree”), the expected PD loss approaches the fraction of extinct species (Fig. 4G). Considering ranked tree shapes shows, however, that the order of nodes has a significant influence on the loss of PD, and on the effect of $$\beta$$ on this loss. If the age and richness of clades are positively correlated ($$\alpha>0$$), the loss of PD is reduced, especially at intermediate extinction fractions (Fig. 4A–F). This is because the smallest subtrees, more prone to early extinction, are younger and hence contain a lower fraction of the PD (Fig. 2). If the age and richness of clades are negatively correlated ($$\alpha<0$$), the loss of PD rises, especially at intermediate extinction fractions. The smallest subtrees, prone to extinction, are older and hence contain more evolutionary distinct species (Fig. 2). This generates losses of PD similar to those observed when the tree shapes are very unbalanced (PD loss equal to the fraction of extinct species). As expected, the effect of $$\alpha$$ is evened out in very unbalanced trees ($$\beta$$ close to $$-2$$; Fig. 4G,H), for which the loss of PD remains close to its highest value whatever the value of $$\alpha$$. In the case of the maximally unbalanced tree shape, there is only one ranked tree shape and the order of nodes is fixed. All these effects of ranked tree shapes on the loss of PD are qualitatively conserved if node depths are distributed as in the Kingman coalescent (instead of the Yule process). In the case of Yule trees, PD loss slightly increases with the initial size of the tree, an effect which is due to more efficient sampling of large values in the common (exponential) distribution of node depths. Yet the results presented above are qualitatively conserved if the size of phylogenetic trees changes (analyses performed with number of species $$N=50$$ and $$N=200$$; see the Supplementary Appendices 3 and 4 available on Dryad). Influence of Nonrandom Extinction Risks on PD Loss The strength of the relation parameterized by $$\eta$$, between the richness of a clade and its relative abundance (here directly influencing the extinction risk of its species) may have a paramount influence on the loss of PD in the face of extinctions (Fig. 6). In trees with ranked tree shapes similar to Yule trees ($$\beta=0$$, $$\alpha=1$$), the concentration of high extinction risks in small clades ($$\eta>1$$) increases the loss of PD, by promoting the extinction of entire clades (Fig. 5). In contrast, when extinction risks are higher in species-richer clades ($$\eta<1$$), phylogenetic redundancy (and hence the likelihood of conserving at least one species per subtree) limits the loss of PD until high extinction levels. Figure 5. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss in Yule trees, for increasing fractions of species extinctions $$p$$. Results are shown either a) as a function of the extinction fraction $$p$$ (for different $$\eta$$ values, with dotted lines showing the bisector) or b) as a function of $$\eta$$ (for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Parameter values: $$\beta=0$$, $$\alpha=0$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 5. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss in Yule trees, for increasing fractions of species extinctions $$p$$. Results are shown either a) as a function of the extinction fraction $$p$$ (for different $$\eta$$ values, with dotted lines showing the bisector) or b) as a function of $$\eta$$ (for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Parameter values: $$\beta=0$$, $$\alpha=0$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 6. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss, for different ranked tree shapes and increasing fractions of species extinctions. Tree balance $$\beta$$ ranges from 10 (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”), and age-richness index $$\alpha$$ ranges from $$-2$$ (A, C, E, G) to $$2$$ (B, D, F, H). Results are shown either as a function of the extinction fraction $$p$$ (left side; for different $$\eta$$ values, and with dotted lines showing the bisector) or as a function of $$\eta$$ (right side; for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 6. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss, for different ranked tree shapes and increasing fractions of species extinctions. Tree balance $$\beta$$ ranges from 10 (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”), and age-richness index $$\alpha$$ ranges from $$-2$$ (A, C, E, G) to $$2$$ (B, D, F, H). Results are shown either as a function of the extinction fraction $$p$$ (left side; for different $$\eta$$ values, and with dotted lines showing the bisector) or as a function of $$\eta$$ (right side; for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. The effect of $$\eta$$ is modified by the ranked shape of species trees. The strength of the age-richness relation (set by $$\alpha$$) modulates the additional loss of PD induced by $$\eta>1$$ (i.e., lower abundances in smaller clades; Fig. 6A–F). When $$\alpha<0$$, smaller clades are not only more prone to extinction but also have deeper nodes, hence more evolutionary distinct species, which increases even further the loss of PD. Unlike in the field of bullets model, the expected PD loss as a function of the fraction $$p$$ of extinct species can even change from convex to concave, and so take values larger than $$p$$ (Fig. 6C,E). When $$\alpha>0$$, smaller clades are more prone to extinction but have shallower nodes, which counteracts the increase of PD loss due to $$\eta>1$$. To summarize, PD loss is increased when $$\eta>1$$ compared to $$\eta=1$$, with a maximal effect for negative values of $$\alpha$$, progressively flattening as $$\alpha$$ grows. We call the “danger zone” the region of parameters corresponding to the theoretical phylogenies that suffer a maximal rate of PD loss straight from the first few extinction events, that is, close to 1% of PD lost for the first 1% of species lost. In the plane $$(\alpha, \eta)$$, the “danger zone” corresponds to $$\{\alpha<0, \eta >1\}$$. As testified by Fig. 6, phylogenies in this zone can even suffer a rate of PD loss which is larger than 1:1 from the first extinction and sustains itself above 1:1 throughout the extinction crisis. In contrast, $$\alpha$$ has little effect on the decrease in PD loss induced by $$\eta<1$$ (i.e., higher abundances in small clades). Indeed, when $$\eta<1$$, the deepest nodes are always protected regardless of the value of $$\alpha$$: when $$\alpha<0$$ the deepest nodes are in small clades which are protected from extinctions by the high relative abundances of its species (due to $$\eta<1$$); when $$\alpha>0$$, the deepest nodes are in large clades which are protected by phylogenetic redundancy. The influence of $$\eta$$ on PD loss is amplified by unbalanced tree shapes ($$\beta<0$$; Fig. 6E–H) and buffered by balanced tree shapes ($$\beta>0$$; Fig. 6A,B), because lower values of $$\beta$$ enhance richness inequalities between clades and raise in turn the influence of $$\eta$$ on PD loss. This interaction between parameters $$\eta$$ and $$\beta$$ overwhelms the influence of $$\alpha$$ (Fig. 6). In the plane $$(\beta, \eta)$$, the “danger zone” is $$\{\beta<-1, \eta >1\}$$ and the previous remark thus implies that in the 3D parameter space, the danger zone is $$\{\alpha <0$$ or $$\beta<-1; \eta >1\}$$. Interestingly, the effect of $$\beta$$ is highly dependent on how extinction risks are distributed within the phylogeny (Fig. 7, and results with other $$\alpha$$ values in Supplementary Appendix S7 available on Dryad). For $$\eta=1$$, we recover the well-known pattern of decreased PD loss as the tree gets more balanced. However, for $$\eta<1$$ we see the reverse pattern, that is PD loss increases with the balance of the tree. Recall that $$\eta <1$$ buffers PD loss, because extinction risks are clustered in the species-richer clades which also display higher phylogenetic redundancy (smaller pendant edges). When the tree is maximally unbalanced, $$\eta<1$$ causes the longest pendant edge to subtend the tip with the largest abundance (and hence to be the last to become extinct). Therefore, the order of extinctions coincides exactly with the increasing order of pendant edge lengths, which results in minimal PD loss for any given level of extinction. In a more balanced phylogeny, the distribution of clade sizes is more even and the buffering effect of the clustered extinction on PD loss is reduced. Figure 7. View largeDownload slide Effect of tree balance $$\beta$$ on PD loss, for different values of clades abundance-richness indexes $$\eta$$, and increasing fractions of species extinctions $$p$$. The clades abundance-richness index $$\eta$$ ranges from $$0.2$$ (left) to 3 (right), and the extinction fraction $$p$$ increases from 0.01 to 0.98 (from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: clades age-richness index $$\alpha=-2$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 7. View largeDownload slide Effect of tree balance $$\beta$$ on PD loss, for different values of clades abundance-richness indexes $$\eta$$, and increasing fractions of species extinctions $$p$$. The clades abundance-richness index $$\eta$$ ranges from $$0.2$$ (left) to 3 (right), and the extinction fraction $$p$$ increases from 0.01 to 0.98 (from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: clades age-richness index $$\alpha=-2$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. For $$\eta >1,$$ we again recover the well-known pattern of decreased PD loss with increasing $$\beta$$. However, when we also have $$\alpha<0$$, the relationship between PD loss and $$\beta$$ is not monotonic, that is for any particular level of extinction, the maximal PD loss is reached for trees with intermediate balance. Recall that $$\alpha<0$$ causes small clades to be relatively older and so to contribute more to PD. The maximal loss of PD thus occurs when extinction risks cluster in small clades. And indeed, when $$\eta>1$$, at each splitting event the species-richer subtree gets a bigger abundance than the species-poorer subtree. However, within a given clade, the abundance of a species should decrease with the number of nodes (splitting events) on its lineage. This latter effect is stronger in unbalanced trees; in balanced trees, extinction risks cannot cluster in small clades, due to the absence of small clades. Trees with intermediate balance do display small clades, and these small clades are large enough to share their low abundance ($$\eta >1$$) into a few species with very low abundance. These species go extinct first, resulting in maximal PD loss. Effect of Species Extinctions on Tree Shape We study the effect of species extinctions on tree shape, seeking in particular to check if the influence of $$\eta$$ on the patterns of PD loss can be explained by changes in tree shape as species become extinct. Figure 8 shows the balance (defined here as the maximum likelihood estimate $$\hat\beta$$ of the parameter $$\beta$$) of the species tree estimated after a fraction $$p$$ of its species have become extinct. When $$\eta=1$$, tree balance is very little altered by extinctions except in very balanced trees, as predicted by the sampling consistency of the model ($$\eta=1$$ amounts to removing species at random except when $$\beta\gg 1$$, see Supplementary Appendix S1 available on Dryad). When $$\eta<1$$, trees tend to become more and more balanced as $$p$$ increases ($$\hat\beta$$ increases with $$p$$), whereas when $$\eta>1$$ trees tend to become more and more similar to Yule trees ($$\hat\beta\to 0$$ as $$p\to 1$$). The effect of $$\eta$$ on PD loss cannot be reduced to its effect on changes in tree shape due to extinctions. On the one hand, $$\eta$$ mostly affects the shape of trees with $$\beta>-1$$ (Fig. 8), whereas tree shape has most effect on PD loss when $$\beta$$ varies between $$-2$$ and $$-1$$ (Fig. 4A,C,E with $$\alpha=0$$). In addition, if the effect of $$\eta$$ on tree shape had a significant influence on PD loss, $$\eta>$$1 should increase this loss when $$\beta>0$$ (by decreasing the balance of trees; Fig. 8D) and decrease it when $$\beta<0$$ (by increasing the balance of trees). Yet, the changes we observe in the effect of $$\eta>1$$ on PD loss for different $$\beta$$ values are the reverse of this prediction. Therefore, the indirect effects of $$\eta$$ (through changes in tree shape) are negligible compared to its direct effects (through nonrandom distribution of extinction risks). Figure 8. View largeDownload slide Effect of abundance-richness index $$\eta$$ on the balance of phylogenetic trees after extinctions (Maximum Likelihood Estimate $$\hat\beta$$ of $$\beta$$). Initial tree balance $$\beta$$ ranges from 10 (brown dots and lines, “bush trees”) to $$-1.9$$ (green dots and lines, “caterpilar trees”). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$, $$\alpha=0$$. Figure 8. View largeDownload slide Effect of abundance-richness index $$\eta$$ on the balance of phylogenetic trees after extinctions (Maximum Likelihood Estimate $$\hat\beta$$ of $$\beta$$). Initial tree balance $$\beta$$ ranges from 10 (brown dots and lines, “bush trees”) to $$-1.9$$ (green dots and lines, “caterpilar trees”). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$, $$\alpha=0$$. As reported above results on the effects of nonrandom extinctions on the loss of PD are conserved when node depths are distributed as in the Kingman coalescent, or when the size of phylogenetic trees changes (analyses performed with $$N=50$$ and $$N=200$$; see Supplementary Appendices S5 and S6 available on Dryad). Parameter Inference When tested against simulated data, the Monte-Carlo inference algorithm by data augmentation performs reasonably well on phylogenies with more than 50 tips for a wide range of parameters (see Supplementary Appendix S8 available on Dryad). As expected, the estimation of $$\beta$$ on trees with at least 50 tips is accurate, since the likelihood formula of the unranked tree is explicit, and this accuracy increases as $$\beta$$ decreases. The inference algorithm also returns overall good estimates of $$\eta$$ and $$\alpha$$ whenever $$\eta>0.3$$. The inference of $$\alpha$$ is unbiased except in the cases where $$\beta < 0$$ and $$\eta < 0.3$$. This corresponds to cases where the unsampled clades are numerous because $$\beta$$ is small, and they have a strong impact on the reconstruction of intervals because $$\eta$$ is small. The inferred $$\eta$$ is overestimated for trees with only $$50$$ tips. For $$\beta < 0$$ and $$\alpha \ge 0$$, $$\eta$$ is slightly overestimated whatever the tip number. For $$\beta > 0$$ and $$\alpha \leq 0,$$ inferences are good for trees with at least $$100$$ tips. Empirical Values Estimates of parameter values on real data shows consistent patterns across all bird family trees. Unsurprisingly, we find negative $$\beta$$ values, mostly comprised between 0 and $$-1$$, corresponding to unbalanced trees (see Supplementary Appendix S9 available on Dryad). Since the estimation of $$\beta$$ is quite accurate for low true values of $$\beta$$ and is biased towards larger estimates than the true value otherwise, these estimates can be taken with confidence. The estimates of $$\eta$$ vary between $$1$$ and $$1.5$$. This indicates that, within bird families, species in small clades tend to have smaller range sizes than species in larger clades. The above study showed that low $$\eta$$ values can be difficult to detect in unbalanced trees. Yet when this is the case, $$\eta$$ is found to be close to the maximal value allowed in the inference (here $$10$$), which is not the case here. We can therefore be confident that these values do not reflect a bias in the inference, but reflect a true pattern in the distribution of range sizes within the phylogenies. Finally, the estimates of $$\alpha$$ are clustered around $$0$$, indicating that there is no correlation between clade sizes and clade depths within each bird phylogeny. This in contrast with what is expected in most explicit models of diversification, where larger clades take more time to diversify, resulting in a strong positive correlation between the depth and the size of clades. When jointly inferring $$\alpha$$ and $$\eta$$ the choice to use range size to infer $$\eta$$ is likely to have an impact on the inferred $$\alpha$$ (because the values of the intervals are reconstructed using tip values, inappropriate tip values would lead to incorrect $$\alpha$$). Therefore, we also ran the inference of $$\alpha$$: we find fairly similar results between values obtained with the inference of $$\alpha$$ only compared to the full inference (the median of the inferred $$\alpha$$ for trees with at least $$50$$ tips is $$0.19$$ when $$\alpha$$ is inferred alone and $$0.05$$ when both $$\alpha$$ and $$\eta$$ are inferred), indicating that tree shape is indeed driving the result (Supplementary Appendix S9 available on Dryad equivalent as Fig. 9 with the $$\alpha$$ inferred without knowledge of the tip range sizes). Figure 9. View largeDownload slide Inferred model parameters on bird family trees of 50 tips or more. Maximum posterior estimate of age-richness index $$\alpha$$ (x-axis), maximum posterior estimate of abundance-richness index $$\eta$$ (y-axis) and maximum likelihood estimate of tree balance $$\beta$$ (point color). Point sizes are proportional to the number of tips in trees, $$N$$. The dashed vertical line shows the value of $$\alpha$$ for trees generated by a birth–death model, and the dashed horizontal line shows the value of $$\eta$$ for which extinction probabilities are distributed within the tree as in a field of bullets model. For all inferences, approximation parameter $$\varepsilon$$ was set to $$0.001$$. Figure 9. View largeDownload slide Inferred model parameters on bird family trees of 50 tips or more. Maximum posterior estimate of age-richness index $$\alpha$$ (x-axis), maximum posterior estimate of abundance-richness index $$\eta$$ (y-axis) and maximum likelihood estimate of tree balance $$\beta$$ (point color). Point sizes are proportional to the number of tips in trees, $$N$$. The dashed vertical line shows the value of $$\alpha$$ for trees generated by a birth–death model, and the dashed horizontal line shows the value of $$\eta$$ for which extinction probabilities are distributed within the tree as in a field of bullets model. For all inferences, approximation parameter $$\varepsilon$$ was set to $$0.001$$. Discussion A New Integrative Measure of the Age-Richness Relation, $$\alpha$$ We introduced here a new model for random ranked tree shapes with a fixed, arbitrary number of tips. This model features two parameters, $$\beta$$ and $$\alpha$$ tuning respectively the balance of the tree and the relative depths or ages of its nodes. Trees with $$\beta\le 0$$ are unbalanced and trees with $$\beta>0$$ are balanced. Whatever the value of $$\alpha$$, the shape of the tree is the same as in Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001). Large clades coalesce deep in the tree when $$\alpha>0$$ and are shallower than smaller clades when $$\alpha<0$$. When $$\beta=0$$ and $$\alpha=1$$, the tree has the same ranked shape as the Kingman coalescent and the Yule tree. Our model does not explicitly incorporate biological processes but enables to generate a broad range of tree shapes by decoupling tree balance and the relation between age and size of clades. In addition, this model is the first model (except the two aforementioned models and the trivial case of the “caterpilar tree”) for ranked tree shapes satisfying sampling-consistency, in the sense that a tree with $$n$$ tips has the same distribution as a tree with $$n+1$$ tips with one tip removed at random. This property is essential to ensure the robustness of the model with respect to incomplete taxon sampling (Heath et al., 2008; Cusimano et al., 2012; Stadler, 2013). Predictions from this model highlight the importance of accounting for node ranks to understand forthcoming changes in macroevolutionary patterns of PD. They show in particular that the relationship between the species richness of a clade and its relative depth in the tree, set by parameter $$\alpha$$ in the model, can have profound impacts on the rate of PD loss (Fig. 4). This parameter $$\alpha$$ constitutes a new index quantifying the age-richness relation of subclades within a phylogeny. A large number of studies have already considered the age-size correlation, assessing its existence (significance, sign, and pattern) across multiple phylogenetic trees, based on one value of species richness and crown or stem age per phylogeny (e.g., Magallon and Sanderson, 2001; Bokma, 2003; Ricklefs, 2006; McPeek and Brown, 2007; Rabosky et al., 2012; Sánchez-Reyes et al., 2016). These studies notably aimed at testing the hypothesis of time-limited diversity patterns, versus hypotheses of diversity set by diversification rates or by limits to diversity (McPeek and Brown, 2007; Ricklefs, 2007; Rabosky, 2009; Barraclough, 2010; Sánchez-Reyes et al., 2016). Our new index $$\alpha$$ is different in that it can be measured by maximizing the likelihood on a single phylogeny, implicitly integrating over all subclades of this phylogeny. An interesting consequence is that one does not have to choose which clades to include in the analysis. For example, $$\alpha$$ is not sensitive to the definition of higher taxa (Stadler et al., 2014; Sánchez-Reyes et al., 2016). Moreover, similarly to what can be done with the index $$\beta$$ (compared to other measures of tree balance; Kirkpatrick and Slatkin, 1993; Aldous, 1996, 2001), a measure of the age-richness relation in a given phylogeny is provided by the maximum likelihood estimate of the model-based parameter $$\alpha$$. Last, we stress that our model does not require the precise knowledge of node datings in the phylogeny but only the relative positions of nodes in time, which preserves $$\alpha$$ estimates from the inaccuracies of time calibrations (Kumar, 2005; Welch and Bromham, 2005; Pulquério and Nichols, 2007; Forest, 2009; Schwartz and Mueller, 2010). Ranked Tree Shapes and the Loss of PD Our results confirm that in the field of bullets model very unbalanced trees undergo stronger loss of PD than balanced trees, under equal fraction of species extinctions. This property was already well known (Nee and May, 1997) but is important to recall given the predominance of unbalanced phylogenetic trees in nature ($$\beta$$ values being often close to $$-1$$; e.g., Guyer and Slowinski, 1991; Heard, 1992; Guyer and Slowinski, 1993; Slowinski and Guyer, 1993; Mooers, 1995; Purvis, 1996; Mooers and Heard, 1997; Blum and François, 2006). However, our results also show that the temporal order of nodes among subtrees (set by the parameter $$\alpha$$) may have even stronger effects than tree balance (set by parameter $$\beta$$; compare the effect of the latter in Supplementary Appendix S8 available on Dryad to that of $$\alpha$$ on Fig. 4). Besides, $$\alpha$$ values below $$0$$ cause drops of PD almost as abrupt as those observed with “caterpilar” shapes ($$\beta$$ close to $$-2$$, with $$\alpha=1$$; Fig. 4D,H). It is therefore essential to consider the ranked shapes of species trees to understand the expected patterns of loss of PD. Values of $$\alpha$$ deviating from $$1$$ may arise from differences in stages of diversification among subtrees, resulting from heterogeneity in biotic or abiotic factors acting on diversification processes in different parts of the species tree. This could be due to bursts of diversification in certain subtrees (e.g., following from key innovations or from migration to empty spatial or ecological space), either recently (resulting in $$\alpha<0$$) or early in the history of clades (resulting in $$\alpha>0$$). Alternatively, $$\alpha$$ values deviating from $$1$$ could be linked to changes in extinction rates in distinct parts of the tree (e.g., due to changes in the biotic or abiotic environment of phylogenetically related species sharing similar ecological niches). Age-dependent speciation (Venditti et al., 2010; Hagen et al., 2015) and extinction (Pearson, 1992; Alexander et al., 2015) are also likely to make node ranking deviate from what is expected in a homogeneous birth–death model. Heterogeneity in diversification rates across the species tree associated with asymmetric competition among species (e.g., evolutionary advantage to previously established species) could limit diversification in younger subtrees, hence leading to $$\alpha>0$$. Last, $$\alpha$$ can be found negative due to the presence of relictual lineages, that is, old clades harboring few species surviving to the present. Modeling Nonrandom Extinctions: $$\eta$$ and the Loss of PD The incorporation of the abundance-richness index $$\eta$$ within the framework provided by Aldous’ $$\beta$$-splitting model allowed us to go beyond the field of bullets assumption. The model allows one to simulate a trait covarying with the shape of the phylogeny. We have so far interpreted this trait as species relative abundance or relative range size, but it could be any species trait with values summing to $$1$$ at the level of the phylogeny. In particular, the trait value of a species can be interpreted as a probability of sampling this species, which is consistent with the initial interpretation of the trait as a relative abundance. In passing, we devised a model of abundance distributions (equivalently interpreted as range size distributions) covarying with the phylogeny, in the broken-stick tradition (MacArthur, 1957; MacArthur and Wilson, 1967). When $$\eta>1$$, the most abundant species are in species-rich clades whereas when $$\eta<1$$ the most abundant species are in species-poor clades. When $$\eta=1$$ all species have the same abundance on average. Here, extinctions are assumed to occur sequentially in the order of increasing abundances. By doing this, we only consider the part of extinction risks that is due to the rarity of a species. In nature, relative extinction risk indeed depends on species abundance, but also on many other features (e.g., dynamics of population growth or decline, fragmentation into subpopulations, biotic or abiotic changes; IUCN, 2012), and may have a significant stochastic component. The simple framework we use to determine extinctions allows us to focus on the direct impact of the distribution of ranked abundances within trees on the loss of PD. This framework can easily be modified to include extrinsic causes of extinctions by taking extinction risk to be a function of abundance and additional known factors. Previous studies concluded that PD loss is increased if extinction risks are clustered in the phylogeny (Davies and Yessoufou, 2013), but that this effect is not substantial (Parhar and Mooers, 2011). Our model shows that the effect on PD loss depends on the way these extinction risks are distributed among clades: PD loss is increased by $$\eta>1$$ (i.e., higher extinction risks in small clades; Fig. 6). Such a distribution of extinction risks may arise from subtrees having low species richness because of higher extinction rates, either due to intrinsic factors (species features that would make them more susceptible to extinction; e.g., long generation time, or low variance or phenotypic plasticity of key ecological traits providing resistance to perturbations or evolutionary advantages in relation to biotic interactions; Purvis et al., 2000c; Johnson et al., 2002) or to extrinsic factors (threats affecting the spatial or ecological space shared by species of the subtree; e.g., Russell et al., 1998; Hughes, 1999; Purvis et al., 2000c; von Euler, 2001; Johnson et al., 2002). Higher extinction risks in small subtrees could also be due to resource limitation affecting simultaneously the density of individuals and the diversity of species, and hence demographic stochasticity (the island effect, Mooers et al., 2009); or to stabilizing selection (e.g., due to competition or to the absence of available spatial or ecological space in the surrounding environment), limiting adaptation and increasing species vulnerability in the face of perturbations (Purvis et al., 2000c; Purvis, 2008). In contrast, $$\eta<1$$ buffers the loss of PD. Higher extinction risks in larger subtrees could result from a trade-off between species richness and average species abundance, provided constrained metacommunity size (with variation along this trade-off following for instance from landscape structure and dynamics, such as geographical isolation affecting the occurrence of allopatric speciation events), from recent speciation events associated with a decrease in average species abundance, geographical range or niche width, or from recent extinction events that removed the most extinction-prone species from certain clades (leaving the latter smaller and with less extinction-prone species; Schwartz and Simberloff, 2001; Lozano and Schwartz, 2005). Hence, $$\eta$$ is expected to vary across clades according to the metacommunity structure and the underpinning diversification dynamics. Given its striking effects on PD loss, this factor should also be accounted for to understand potential future losses of PD. Combined Effects of $$\beta$$, $$\alpha$$, and $$\eta$$: Reversing Some Expected Patterns of PD Loss The influence of $$\eta$$ on the loss of PD is enhanced by $$\alpha<0$$ (small clades containing evolutionary distinct species) and $$\beta<0$$ (more variability in clades richness) (Fig. 6). However, a stronger clustering of extinction risks does not necessarily lead to higher loss of PD (e.g., if extinctions occur first in species-richer subtrees—which contain more phylogenetic redundancy—as in the case when $$\alpha>0$$ and $$\eta<1$$). These interactions between the effects of $$\beta$$, $$\alpha,$$ and $$\eta$$ may reverse two well-known patterns of variation in the loss of PD (Nee and May, 1997). First, the increase in PD loss with tree imbalance can be hampered by $$\eta$$ values deviating from one (Fig. 7 and Supplementary Appendix S7 available on Dryad). In particular when $$\eta<1$$, this pattern results from the preferential extinction of phylogenetically redundant species in more unbalanced trees when extinction risks are clustered in large clades. Second, when $$\eta>1$$ and $$\alpha<0$$ the loss of PD proceeds faster than that of species diversity (turning their relationship from convex to concave, except in very balanced or very unbalanced trees; Fig. 6C,E). This pattern (also highlighted in Heard and Mooers, 2000) is caused by the preferential extinction in small subtrees containing evolutionary distinct species. The only other cases where such high loss of PD is reached is when $$\beta<-1$$ and $$\eta>1$$. This parameter zone is a macroevolutionary danger zone, where particular phylogenetic shapes combine with clustered patterns of extinction probabilities to produce large and rapid losses of evolutionary history, much like ice and wind bringing down thick branches of trees in winter. In this region of the parameter space ($$\beta <-1$$ or $$\alpha <0$$; $$\eta>1$$), phylogenies are prone to a sudden collapse of PD. Loss of PD in Bird Family Phylogenies Our inference study shows that the phylogeny of bird families tend to exhibit $$\beta$$ values comprised between $$-1$$ and $$0$$. A similar result was found in many macroevolutionary studies, commonly observing values of $$\beta$$ clustering around $$-1$$ in real phylogenies (e.g., Guyer and Slowinski, 1991; Heard, 1992; Guyer and Slowinski, 1993; Slowinski and Guyer, 1993; Mooers, 1995; Purvis, 1996; Mooers and Heard, 1997; Blum and François, 2006). With these topologies, we expect both $$\alpha$$ and $$\eta$$ to play a major role in determining the potential losses of PD (Fig. 6E,F). We observed $$\alpha$$ values clustering around zero, consistently with several empirical studies that found no positive relation between clade age and clade size (Ricklefs, 2007, 2009; Rabosky et al., 2012). These values contrast with the value of $$1$$ expected in Yule trees, and make phylogenies very sensitive to PD loss. Our estimates of $$\eta$$ values, based on the distribution of range sizes in bird family phylogenies, all fall between $$1$$ and $$1.5$$. This indicates that species in small clades tend to have smaller ranges than species in bigger clades. Range size has been shown to be one of the most important correlates of extinction risks and is one of the IUCN red list criteria (Purvis et al., 2000b; Cardillo et al., 2006; Lee and Jetz, 2011; IUCN, 2012; Arbetman et al., 2017). Expectedly, range sizes in our data set are significantly different among IUCN status classes (see Supplementary Appendix S9 available on Dryad). This pattern of extinction risks clustering in species-poor clades, highlighted in our data by $$\eta$$ values above $$1$$, has also been shown in plants (Vamosi and Wilson, 2008), and in birds and mammals for past extinctions (Russell et al., 1998; Mooers et al., 2009). Considering the three parameters together, we find that bird family trees are situated close to the region of the parameter space termed “danger zone,” for which we find the loss of PD to be at least as fast as the loss of species diversity. In particular, the combination of negative $$\alpha$$ values with $$\eta>1$$ leads to higher extinction risks for evolutionary distinct species. We can expect such a pattern as a result from evolutionary mechanisms acting simultaneously on different features of trees. For example, subtree-specific susceptibility to extinction, or stabilizing selection generating relictual lineages, are both expected to beget small subtrees with high divergence times also endowed with high species extinction risks. This pattern has been already found for past extinctions in birds using a species level measure of evolutionary distinctiveness; authors observed in that case a similar loss of species and PD (von Euler, 2001; Szabo et al., 2012). Evolutionary distinct bird lineages were also shown to be more threatened by agricultural expansion and intensification than more recent lineages in Costa Rica (Frishkoff et al., 2014). This was also found in other taxa, such as marsupial mammals (Johnson et al., 2002) and Sebastes (Magnuson-Ford et al., 2009) (but see Davies et al. 2011, who show that in South African plants, extinction risks cluster in younger, fast-diversifying genera). A striking result of our inference study relates to the narrow range of $$\alpha$$ values obtained as soon as the trees are large enough for the inference to be accurate (see Supplementary Appendix S9 available on Dryad for inferred parameter values as a function of the tip number in the phylogenies). This value, which differs from what is found in birth–death models, adds a new conundrum concerning the shape of empirical trees. Branch Lengths in Empirical Phylogenies The parameter $$\alpha$$ of the model shapes the order in which speciations take place, but does not instantiate the actual times between two consecutive speciation events, that is, edge lengths. In the numerical investigations of PD loss, we considered two models for edge lengths: the pure-birth process (Yule, 1925) and the Kingman coalescent (Kingman, 1982). Using either of these models did not affect our results qualitatively, but affected them quantitatively (compare Figs. 4 and 6 to figures provided in Supplementary Appendices S4 and S6 available on Dryad). Our modeling framework allows easy exploration of predictions under different models of edge lengths. This is interesting as many empirical phylogenies are not time-calibrated, or imprecisely. Besides, empirical phylogenetic trees were shown to often exhibit a decrease in branching tempo, that is, in the rate of lineage accumulation through time (characterized in particular by estimates of the statistic $$\gamma<0$$; e.g., Nee et al., 1992; Zink and Slowinski, 1995; Lovette and Bermingham, 1999; Pybus and Harvey, 2000; Rüber and Zardoya, 2005; Kozak et al., 2006; Seehausen, 2006; Weir, 2006; McPeek, 2008; Phillimore and Price, 2008; Rabosky and Lovette, 2008a; Jønsson et al., 2012; Moen and Morlon, 2014). Hence, quantitative predictions on the loss of phylogenetic diversity in the face of species extinctions could be further increased by accounting for real branch lengths. Moreover, several theoretical studies suggested that the branching tempo of species trees may change with clade age, decreasing in particular in younger clades (the “out of equilibrium” hypothesis, proposed to explain the negative values of $$\gamma$$ often observed in real phylogenies; Liow et al., 2010; Gascuel et al., 2015; Manceau et al., 2015; Missa et al., 2016; Bonnet-Lebrun et al., 2017). Taking into account such correlations between the age of clades and their branching tempo would also affect the expected loss of PD. The EDGE program (“Evolutionary Distinct and Globally Endangered”; Isaac et al., 2007) encourages conservation priorities aiming at preserving most evolutionary history within the tree of life, by proposing a ranking of species based on combined criteria of evolutionary distinctiveness and extinction risk. Although our approach is not species-based but clade-based, it also investigates the preservation of evolutionary history based on principles linked to species evolutionary distinctiveness (related to the depths of subtrees, which depend on $$\alpha$$) and to the distribution of extinction risks in the tree (which depends on $$\eta$$). Accordingly, to conserve most evolutionary history and evolutionary potential for further diversification and/or survival, priority could be given to clades that would undergo higher loss of PD in the face of species extinctions, that is, clades in the danger zone ($$\eta>1$$ and either $$\beta<-1$$ or $$\alpha<0)$$, and although not shown but only discussed herein, with $$\gamma<0$$ (decreasing branching tempo; Pybus and Harvey, 2000). Beyond Losses of PD As we have seen earlier, the trait value of a species initially interpreted as its relative abundance can also be seen as a probability of sampling it. So an additional tool provided by the parameter $$\eta$$ is a sampling distribution where sampling is not independent of the phylogeny. In particular, our results could be interpreted in the light of rarefaction experiments (Nipperess and Matsen, 2013), which study the way phylogenetic patterns in a metacommunity change as sampling decreases. Previous studies already pointed out strong impacts of nonrandom taxon sampling on the macroevolutionary patterns that we observe (e.g., Cusimano and Renner, 2010). Our results provide insights on the effects of nonrandom sampling on PD and phylogenetic tree topology. They reveal how, when the rarest species are not sampled, the discrepancy between observed and real PD depends on tree balance and abundance-richness relation (this discrepancy being larger in particular in the danger zone; Fig. 6E); and how the discrepancy between observed and real tree shape depends on $$\eta$$ (real trees being more unbalanced if $$\eta<1$$, and diverging from Yule trees towards more balance or more unbalance if $$\eta>1$$; Fig. 8). These effects of incomplete sampling on macroevolutionary patterns should be particularly important to understand biodiversity patterns in bacterial and archeal phyla, which remain poorly known in particular because they likely harbor rare species having high chances to remain unnoticed. Conclusion This new stochastic model of phylogenetic trees spans a large range of binary trees endowed with node rankings and species abundances/range sizes/extinction risks, based on three parameters only and interpolating other well-known models. We showed that ranked tree shapes, nonrandom extinctions and the interactions thereof, may have a strong impact on the loss of PD in the face of species extinctions, potentially reversing some expected patterns of variation in PD. The simplicity of the model allows one to infer the parameters on empirical phylogenies. Applying our inference procedure on bird family phylogenies we found that, in this data set, the parameters fall within a narrow range of the parameter space; and that the inferred values make the PD of these trees very sensitive to species extinctions. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.mv980. Acknowledgments The authors are very grateful to Mike Steel, Ana S.L. Rodrigues, and Hélène Morlon for their feedbacks on an earlier version of this manuscript. They are indebted to Arne Mooers and to Walter Jetz for providing the data on bird range sizes. They wish to thank the Associate Editor and the referees for their careful reviews. They thank the Center for Interdisciplinary Research in Biology (Collège de France, CNRS) for funding. References Aldous D. 1996. Probability distributions on cladograms. In: Aldous, D., Pemantle, R., editors. Random discrete structures . New York: Springer. p. 1– 18. Google Scholar CrossRef Search ADS   Aldous D. 2001. Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Stat. Sci.  16: 23– 34. Google Scholar CrossRef Search ADS   Alexander H.K., Lambert A., Stadler T. 2015. Quantifying age-dependent extinction from species phylogenies. Syst. Biol.  65: 35– 50. Google Scholar CrossRef Search ADS PubMed  Alfaro M., Santini F., Brock C., Alamillo H., Dornburg A., Rabosky D., Carnevale G., Harmon L. 2009. Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates. Proc. Natl. Acad. Sci. USA  106: 13410– 13414. Google Scholar CrossRef Search ADS   Arbetman, M.P., Gleiser G., Morales C.L., Williams P., Aizen M.A. 2017. Global decline of bumblebees is phylogenetically structured and inversely related to species range size and pathogen incidence. Proc. Roy. Soc. Lond. B  284: 20170204. Google Scholar CrossRef Search ADS   Baillie, J.E.M., Hilton-Taylor C., Stuart S.N. 2004. A global species assessment. Gland, Switzerland: IUCN. Google Scholar CrossRef Search ADS   Barnosky, A.D., Matzke N., Tomiya S., Wogan G.O.U., Swartz B., Quental T.B., Marshall C., McGuire J.L., Lindsey E.L., Maguire K.C., Mersey B., Ferrer E.A. 2011. Has the Earth’s sixth mass extinction already arrived? Nature  471: 51– 57. Google Scholar CrossRef Search ADS PubMed  Barraclough, T.G. 2010. Evolving entities: towards a unified framework for understanding diversity at the species and higher levels. Philos. Trans. R. Soc. Lond. B  365: 1801– 1813. Google Scholar CrossRef Search ADS   Bennett, P.M., Owens I.P.F. 1997. Variation in extinction risk among birds: chance or evolutionary predisposition? Proc. R. Soc. Lond. B.  264: 401– 408. Google Scholar CrossRef Search ADS   Bertoin, J. 2002. Self-similar fragmentations. Ann. I.H. Poincaré  38: 319– 340. Google Scholar CrossRef Search ADS   Bertoin, J. 2006. Random fragmentation and coagulation processes. Cambridge: Cambridge University Press. Google Scholar CrossRef Search ADS   Bielby, J., Cunningham A.A., Purvis A. 2006. Taxonomic selectivity in amphibians: ignorance, geography or biology? Anim.  Conserv.  9: 135– 143. Blum, M., François O. 2006. Which random processes describe the tree of life? A large-scale study of phylogenetic tree imbalance. Syst. Biol.  55: 685– 691. Bokma, F. 2003. Testing for equal rates of cladogenesis in diverse taxa. Evolution  57: 2469– 2474. Google Scholar CrossRef Search ADS PubMed  Bonnet-Lebrun, A.-S., Manica A., Eriksson A., Rodrigues A.S. 2017. Empirical phylogenies and species abundance distributions are consistent with preequilibrium dynamics of neutral community models with gene flow. Evolution  71: 1149– 1163. Google Scholar CrossRef Search ADS PubMed  Bortolussi N., Durand E., Blum M., François O. 2006. apTreeshape: statistical analysis of phylogenetic tree shape. Bioinformatics  22: 363– 364. Google Scholar CrossRef Search ADS PubMed  Cardillo, M., Mace G.M., Gittleman J.L., Purvis A. 2006. Latent extinction risk and the future battlegrounds of mammal conservation. Proc. Roy. Soc. Lond. B  103: 4157– 4161. Cardillo M., Orme C.D.L., Owens I.P.F. 2005. Testing for latitudinal bias in diversification rates : an example using New World birds. Ecology  86: 2278– 2287. Google Scholar CrossRef Search ADS   Colwell R.K., Lees D.C. 2000. The mid-domain effect: geometric species richness. Trends Ecol. Evol.  15: 70– 76. Google Scholar CrossRef Search ADS PubMed  Cusimano N., Renner S.S. 2010. Slowdowns in diversification rates from real phylogenies may not be real. Syst. Biol.  59: 458– 464. Google Scholar CrossRef Search ADS PubMed  Cusimano N., Stadler T., Renner S.S. 2012. A new method for handling missing species in diversification analysis applicable to randomly or nonrandomly sampled phylogenies. Syst. Biol.  61: 785– 792. Google Scholar CrossRef Search ADS PubMed  Davies T.J., Smith G.F., Bellstedt D.U., Boatwright J.S., Bytebier B., Cowling R.M., Forest F., Harmon L.J., Muasya A.M., Schrire B.D., Steenkamp Y., van der Bank M., Savolainen V. 2011. Extinction risk and diversification are linked in a plant biodiversity hotspot. PLoS Biol.  9: 1– 9. Google Scholar CrossRef Search ADS   Davies T.J., Yessoufou K. 2013. Revisiting the impacts of non-random extinction on the tree-of-life. Biol. Lett.  9: 20130343. Google Scholar CrossRef Search ADS PubMed  Doran N.A., Arnold A.J., Parker W.C., Huffer F.W. 2006. Is extinction age dependent? Palaios  21: 571– 579. Google Scholar CrossRef Search ADS   Ewens W.J. 2012. Mathematical population genetics 1: theoretical introduction. New York: Springer Science & Business Media. Google Scholar CrossRef Search ADS   Faith D.P. 1992. Conservation evaluation and phylogenetic diversity. Biol. Conserv.  61: 1– 10. Google Scholar CrossRef Search ADS   Faller B., Pardi F., Steel M. 2008. Distribution of phylogenetic diversity under random extinction. J. Theor. Biol.  251: 286– 296. Google Scholar CrossRef Search ADS PubMed  Feng S. 2010. The Poisson–Dirichlet distribution and related topics: models and asymptotic behaviors. Heidelberg: Springer Science & Business Media. Google Scholar CrossRef Search ADS   Forest F. 2009. Calibrating the tree of life: fossils, molecules and evolutionary timescales. Ann. Bot.  104: 789– 794. Google Scholar CrossRef Search ADS PubMed  Frishkoff L.O., Karp D.S., M’Gonigle L.K., Mendenhall C.D., Zook J., Kremen C., Hadly E.A., Daily G.C. 2014. Loss of avian phylogenetic diversity in neotropical agricultural systems. Science  345: 1343– 1346. Google Scholar CrossRef Search ADS PubMed  Fritz S.A., Purvis A. 2010. Selectivity in mammalian extinction risk and threat types: a new measure of phylogenetic signal strength in binary traits. Conserv. Biol.  24: 1042– 1051. Google Scholar CrossRef Search ADS PubMed  Gascuel F., Ferrière R., Aguilée R., Lambert A. 2015. How ecology and landscape dynamics shape phylogenetic trees. Syst. Biol.  64: 590– 607. Google Scholar CrossRef Search ADS PubMed  Gaston K.J., Blackburn T.M. 1995. Birds, body size and the threat of extinction. Philos. Trans. R. Soc. Lond. B  347: 205– 212. Google Scholar CrossRef Search ADS   Glavin T. 2007. The sixth extinction: journeys among the lost and left behind. New York: Thomas Dunne Books. Guyer C., Slowinski J.B. 1991. Comparisons of observed phylogenetic topologies with null expectations among three monophyletic lineages. Evolution  45: 340– 350. Google Scholar CrossRef Search ADS PubMed  Guyer C., Slowinski J.B. 1993. Adaptive radiation and the topology of large phylogenies. Evolution  47: 253– 263. Google Scholar CrossRef Search ADS PubMed  Hagen O., Hartmann K., Steel M., Stadler T. 2015. Age-dependent speciation can explain the shape of empirical phylogenies. Syst. Biol.  64: 432– 440. Google Scholar CrossRef Search ADS PubMed  Heard S.B. 1992. Patterns in tree balance among cladistic, phenetic, and randomly generated phylogenetic trees. Evolution  46: 1818– 1826. Google Scholar CrossRef Search ADS PubMed  Heard S.B., Mooers A.O. 2000. Phylogenetically patterned speciation rates and extinction risks change the loss of evolutionary history during extinctions. Proc. R. Soc. Lond. B.  267: 613– 620. Google Scholar CrossRef Search ADS   Heath T.A., Hedtke S.M., Hillis D.M. 2008. Taxon sampling and the accuracy of phylogenetic analyses. J. Syst. Evol.  46: 239– 257. Hubbell S. 2001. The unified neutral theory of biodiversity and biogeography. Princeton, NJ: Princeton University Press. Google Scholar CrossRef Search ADS   Hughes A.L. 1999. Differential human impact on the survival of genetically distinct avian lineages. Bird Conserv. Int.  9: 147– 154. Google Scholar CrossRef Search ADS   Isaac N.J.B., Turvey S.T., Collen B., Waterman C., Baillie J.E.M 2007. Mammals on the EDGE: conservation priorities based on threat and phylogeny. PLoS One  2: e296. Google Scholar CrossRef Search ADS PubMed  IUCN. 2012. IUCN red list categories and criteria: version 3.1. 2nd ed. Gland, Switzerland: IUCN. Jetz W., Thomas G., Joy J., Hartmann K., Mooers A. 2012. The global diversity of birds in space and time. Nature  491: 444– 448. Google Scholar CrossRef Search ADS PubMed  Johnson C.N., Delean S., Balmford A. 2002. Phylogeny and the selectivity of extinction in Australian marsupials. Anim. Conserv.  5: 135– 142. Google Scholar CrossRef Search ADS   Johnson S.G., Narasimhan B. 2013. R package “cubature”: adaptive multivariate integration over hypercubes. Available from: https://cran.r-project.org/web/packages/cubature/index.html. Jønsson, K.A., Fabre P.-h., Fritz S.A., Etienne R.S., Ricklefs R.E., Jørgensen T.B. 2012. Ecological and evolutionary determinants for the adaptive radiation of the Madagascan vangas. Proc. Natl. Acad. Sci. USA  109: 6620– 6625. Google Scholar CrossRef Search ADS   Kembel S.W., Ackerly D.D., Blomberg S.P., Cornwell W.K., Cowan P.D., Helmus M.R., Morlon H., Webb C.O. 2014. R package “picante”: R tools for integrating phylogenies and ecology. Available from: https://cran.r-project.org/web/packages/picante/index.html. Kingman J.F.C. 1982. The coalescent. Stoch. Proc. Appl.  13: 235– 248. Google Scholar CrossRef Search ADS   Kirkpatrick M., Slatkin M. 1993. Searching for evolutionary patterns in the shape of a phylogenetic tree. Evolution  47: 1171– 1181. Google Scholar CrossRef Search ADS PubMed  Kozak K.H., Weisrock D.W., Larson A. 2006. Rapid lineage accumulation in a non-adaptive radiation: phylogenetic analysis of diversification rates in eastern North American woodland salamanders (Plethodontidae: Plethodon). Proc. R. Soc. Lond. B.  273: 539– 546. Google Scholar CrossRef Search ADS   Kumar S. 2005. Molecular clocks: four decades of evolution. Nat. Rev. Genet.  6: 654– 662. Google Scholar CrossRef Search ADS PubMed  Lambert A., Stadler T. 2013. Birth-death models and coalescent point processes: the shape and probability of reconstructed phylogenies. Theor. Popul. Biol.  90: 113– 128. Google Scholar CrossRef Search ADS PubMed  Lambert A., Steel M. 2013. Predicting the loss of phylogenetic diversity under non-stationary diversification models. J. Theor. Biol.  337: 111– 124. Google Scholar CrossRef Search ADS PubMed  Lambert A.et al.   2017. Probabilistic models for the (sub) tree (s) of life. Braz. J. Probab. Stat.  31: 415– 475. Google Scholar CrossRef Search ADS   Leakey R.E., Lewin R. 1995. The sixth extinction: patterns of life and the future of humankind. New York: Doubleday. Lee T.M., Jetz W. 2011. Unravelling the structure of species extinction risk for predictive conservation science. Proc. Roy. Soc. Lond. B.  278 1710: 1329– 1338 Google Scholar CrossRef Search ADS   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: 646– 659. Google Scholar CrossRef Search ADS PubMed  Lovette I.J., Bermingham E. 1999. Explosive speciation in the New World Dendroica warblers. Proc. R. Soc. Lond. B.  266: 1629– 1636. Google Scholar CrossRef Search ADS   Lozano F.D., Schwartz M.W. 2005. Patterns of rarity and taxonomic group size in plants. Biol. Conserv.  126: 146– 154. Google Scholar CrossRef Search ADS   MacArthur R., Wilson E. 1967. The theory of island biogeography. Princeton, NJ: Princeton University Press. Google Scholar CrossRef Search ADS   MacArthur R.H. 1957. On the relative abundance of bird species. Proc. Natl. Acad. Sci. USA  43: 293– 295. Google Scholar CrossRef Search ADS   Maddison W.P., Midford P.E., Otto S.P. 2007. Estimating a binary character’s effect on speciation and extinction. Syst. Biol.  56: 701– 710. Google Scholar CrossRef Search ADS PubMed  Magallon S., Sanderson M.J. 2001. Absolute diversification rates in Angiosperm clades. Evolution  55: 1762– 1780. Google Scholar CrossRef Search ADS PubMed  Magnuson-Ford K., Ingram T., Redding D.W., Mooers A.Ø. 2009. Rockfish (sebastes) that are evolutionarily isolated are also large, morphologically distinctive and vulnerable to overfishing. Biol. Conserv.  142: 1787– 1796. Google Scholar CrossRef Search ADS   Manceau M., Lambert A., Morlon H. 2015. Phylogenies support out-of-equilibrium models of biodiversity. Ecol. Lett.  18: 347– 356. Google Scholar CrossRef Search ADS PubMed  McKinney M.L. 1997. Extinction vulnerability and selectivity: combining ecological and paleontological views. Annu. Rev. Ecol. Syst.  28: 495– 516. Google Scholar CrossRef Search ADS   McPeek M.A. 2008. The ecological dynamics of clade diversification and community assembly. Am. Nat.  172: E270– E284. Google Scholar CrossRef Search ADS PubMed  McPeek M.A., Brown J.M. 2007. Clade age and not diversification rate explains species richness among animal taxa. Am. Nat.  169: E97– 106. Google Scholar CrossRef Search ADS PubMed  Missa O., Dytham C., Morlon H. 2016. Understanding how biodiversity unfolds through time under neutral theory. Philos. Trans. R. Soc. Lond. B  371: 1– 12. Google Scholar CrossRef Search ADS   Moen D., Morlon H. 2014. Why does diversification slow down? Trends Ecol. Evol.  29: 190– 197. Google Scholar CrossRef Search ADS PubMed  Mooers A., Gascuel O., Stadler T., Li H., Steel M. 2012. Branch lengths on birth-death trees and the expected loss of phylogenetic diversity. Syst. Biol.  61: 195– 203. Google Scholar CrossRef Search ADS PubMed  Mooers A., Heard S. 1997. Inferring evolutionary process from phylogenetic tree shape. Q. Rev. Biol.  72: 31– 54. Google Scholar CrossRef Search ADS   Mooers A.O. 1995. Tree balance and tree completeness. Evolution  49: 379– 384. Google Scholar CrossRef Search ADS PubMed  Mooers A. Ø., Goring S.J., Turvey S.T., Kuhn T.S. 2009. Holocene extinctions and the loss of feature diversity. Turvey S. editor. Holocene extinctions.  Oxford: Oxford University Press. p. 279– 338 Google Scholar CrossRef Search ADS   Morlon H., Parsons T.L., Plotkin J.B. 2011a. Reconciling molecular phylogenies with the fossil record. Proc. Natl. Acad. Sci. USA  108: 16327– 16332. Google Scholar CrossRef Search ADS   Morlon H., Schwilk D.W., Bryant J.A., Marquet P.A., Rebelo A.G., Tauss C., Bohannan B.J., Green J.L. 2011b. Spatial patterns of phylogenetic diversity. Ecol. Lett.  14: 141– 149. Google Scholar CrossRef Search ADS   Nee S. 2006. Birth-death models in macroevolution. Annu. Rev. Ecol. Evol. Syst.  37: 1– 17. Google Scholar CrossRef Search ADS   Nee S., May R.M. 1997. Extinction and the loss of evolutionary history. Science  278: 692– 694. Google Scholar CrossRef Search ADS PubMed  Nee S., Mooers A.O., Harvey P.H. 1992. Tempo and mode of evolution revealed from molecular phylogenies. Proc. Natl. Acad. Sci. USA  89: 8322– 8326. Google Scholar CrossRef Search ADS   Nipperess D.A., Matsen F.A. 2013. The mean and variance of phylogenetic diversity under rarefaction. Methods Ecol. Evol.  4: 566– 572. Google Scholar CrossRef Search ADS PubMed  Paradis E., Claude J., Strimmer K. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics  20: 289– 290. Google Scholar CrossRef Search ADS PubMed  Parhar R.K., Mooers A.Ø. 2011. Phylogenetically clustered extinction risks do not substantially prune the tree of life. PLoS One  6: e23528. Google Scholar CrossRef Search ADS PubMed  Pearson P.N. 1992. Survivorship analysis of fossil taxa when real-time extinction rates vary: the paleogene planktonic foraminifera. Paleobiology  18: 115– 131. Google Scholar CrossRef Search ADS   Phillimore A.B., Price T.D. 2008. Density-dependent cladogenesis in birds. PLoS Biol . 6: e71. Google Scholar CrossRef Search ADS PubMed  Pigot A.L., Phillimore A.B., Owens I.P.F., Orme C.D.L. 2010. The shape and temporal dynamics of phylogenetic trees arising from geographic speciation. Syst. Biol.  59: 660– 673. Google Scholar CrossRef Search ADS PubMed  Prado P.I., Miranda M.D., Chalom A. 2015. R package “sads”: maximum likelihood models for species abundance distributions. Available form: http://search.r-project.org/library/sads/html/fitsad.html. Pulquério, M.J.F., Nichols R.A. 2007. Dates from the molecular clock: how wrong can we be? Trends Ecol. Evol.  22: 180– 184. Google Scholar CrossRef Search ADS PubMed  Purvis A. 1996. Using interspecies phylogenies to test macroevolutionary hypotheses. In: Harvey, P., Leigh Brown, A., Maynard Smith, J., Nee, S., editors. New uses for new phylogenies . Oxford: Oxford University Press. p. 153– 168. Purvis A. 2008. Phylogenetic approaches to the study of extinction. Annu. Rev. Ecol. Evol. Syst.  39: 301– 319. Google Scholar CrossRef Search ADS   Purvis A., Agapow P.-M., Gittleman J.L., Mace G.M. 2000a. Nonrandom extinction and the loss of evolutionary history. Science  288: 328– 330. Google Scholar CrossRef Search ADS   Purvis A., Gittleman J.L., Cowlishaw G., Mace G.M. 2000b. Predicting extinction risk in declining species. Proc. Roy. Soc. Lond. B  267: 1947– 1952. Google Scholar CrossRef Search ADS   Purvis A., Jones K.E., Mace G.M. 2000c. Extinction. BioEssays  22: 1123– 1133. Google Scholar CrossRef Search ADS   Pybus O.G., Harvey P.H. 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc. R. Soc. Lond. B.  267: 2267– 2272. Google Scholar CrossRef Search ADS   R Development Core Team. 2012. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Rabosky D.L. 2009. Ecological limits and diversification rate: alternative paradigms to explain the variation in species richness among clades and regions. Ecol. Lett.  12: 735– 743. Google Scholar CrossRef Search ADS PubMed  Rabosky D.L., Lovette I.J. 2008a. Density-dependent diversification in North American wood warblers. Proc. R. Soc. Lond. B.  275: 2363– 2371. Google Scholar CrossRef Search ADS   Rabosky D.L., Lovette I.J. 2008b. Explosive evolutionary radiations: decreasing speciation or increasing extinction through time? Evolution  62: 1866– 1875. Google Scholar CrossRef Search ADS   Rabosky D.L., Slater G.J., Alfaro M.E. 2012. Clade age and species richness are decoupled across the eukaryotic tree of life. PLoS Biol . 10: e1001381. Google Scholar CrossRef Search ADS PubMed  Raup D.M., Gould S.J., Schopf T.J.M., Simberloff D.S. 1973. Stochastic models of phylogeny and the evolution of diversity. J. Geol.  81: 525– 542. Google Scholar CrossRef Search ADS   Redding D.W., Mooers A.Ø. 2006. Incorporating evolutionary measures into conservation prioritization. Conserv. Biol.  20: 1670– 1678. Google Scholar CrossRef Search ADS PubMed  Ricklefs R. 2009. Speciation, extinction and diversity. In: Butlin, R., Bridle, J., Schluter, D., editors. Speciation and patterns of diversity . Cambridge: Cambridge University Press. p. 257– 277. Google Scholar CrossRef Search ADS   Ricklefs R.E. 2006. Global variation in the diversification rate of passerine birds. Ecology  87: 2468– 78. Google Scholar CrossRef Search ADS PubMed  Ricklefs R.E. 2007. History and diversity: explorations at the intersection of ecology and evolution. Am. Nat.  170: S56– S70. Google Scholar CrossRef Search ADS PubMed  Rüber L., Zardoya R. 2005. Rapid cladogenesis in marine fishes revisited. Evolution  59: 1119– 1127. Google Scholar CrossRef Search ADS PubMed  Russell G.J., Brooks T.M., McKinney M.M., Anderson C.G. 1998. Present and future taxonomic selectivity in bird and mammal extinctions. Conserv. Biol.  12: 1365– 1376. Google Scholar CrossRef Search ADS   Sainudiin R., Véber A. 2016. A beta-splitting model for evolutionary trees. R. Soc. Open Sci.  3: 160016. Google Scholar CrossRef Search ADS PubMed  Sánchez-Reyes L.L., Morlon H., Magallón S. 2016. Uncovering higher-taxon diversification dynamics from clade age and species-richness data. Syst. Biol.  66: 367– 378. Schwartz M.W., Simberloff D. 2001. Taxon size predicts rates of rarity in vascular plants. Ecol. Lett.  4 5: 464– 469. Google Scholar CrossRef Search ADS   Schwartz R.S., Mueller R.L. 2010. Branch length estimation and divergence dating: estimates of error in Bayesian and maximum likelihood frameworks. BMC Evol. Biol.  10: 1– 21. Google Scholar CrossRef Search ADS PubMed  Seehausen O. 2006. African cichlid fish: a model system in adaptive radiation research. Proc. R. Soc. Lond. B.  273: 1987– 1998. Google Scholar CrossRef Search ADS   Slowinski J., Guyer C. 1993. Testing whether certain traits have caused amplified diversification—an improved method based on a model of random speciation and extinction. Am. Nat.  142: 1019– 1024. Google Scholar CrossRef Search ADS PubMed  Stadler T. 2013. Recovering speciation and extinction dynamics based on phylogenies. J. Evol. Biol.  26: 1203– 1219. Google Scholar CrossRef Search ADS PubMed  Stadler T., Rabosky D.L., Ricklefs R.E., Bokma F. 2014. On age and species richness of higher taxa. Am. Nat.  184: 447– 455. Google Scholar CrossRef Search ADS PubMed  Szabo J.K., Khwaja N., Garnett S.T., Butchart S.H. 2012. Global patterns and drivers of avian extinctions at the species and subspecies level. PLoS One  7: e47080. Google Scholar CrossRef Search ADS PubMed  Vamosi J.C., Wilson J.R. 2008. Nonrandom extinction leads to elevated loss of angiosperm evolutionary history. Ecol. Lett.  11: 1047– 1053. Google Scholar CrossRef Search ADS PubMed  Van Valen L. 1976. Ecological species, multispecies, and oaks. Taxon  25: 233– 239. Google Scholar CrossRef Search ADS   Vazquez D.P., Gittleman J.L. 1998. Biodiversity conservation: does phylogeny matter? Curr. Biol.  8: 379– 381. Google Scholar CrossRef Search ADS   Venditti C., Meade A., Pagel M. 2010. Phylogenies reveal new interpretation of speciation and the Red Queen. Nature  463: 349– 352. Google Scholar CrossRef Search ADS PubMed  Veron S., Davies T.J., Cadotte M.W., Clergeau P., Pavoine S. 2015. Predicting loss of evolutionary history: where are we? Biol. Rev.  0: 0– 0. von Euler F. 2001. Selective extinction and rapid loss of evolutionary history in the bird fauna. Proc. R. Soc. Lond. B.  268: 127– 130. Google Scholar CrossRef Search ADS   Wake D.B., Vredenburg V.T. 2008. Are we in the midst of the sixth mass extinction? A view from the world of amphibians. Proc. Natl. Acad. Sci. USA  105: 11466– 11473. Google Scholar CrossRef Search ADS   Weir J.T. 2006. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution  60: 842– 855. Google Scholar CrossRef Search ADS PubMed  Welch J.J., Bromham L. 2005. Molecular dating when rates vary. Trends Ecol. Evol.  20: 320– 327. Google Scholar CrossRef Search ADS PubMed  Yule G. 1925. A mathematical theory of evolution, based on the conclusions of Dr. J.C. Willis, F.R.S. Philos. Trans. R. Soc. Lond. B  213: 402– 410. Google Scholar CrossRef Search ADS   Zink R., Slowinski J. 1995. Evidence from molecular systematics for decreased avian diversification in the Pleistocene Epoch. Proc. Natl. Acad. Sci. USA  92: 5832– 5835. Google Scholar CrossRef Search ADS   © 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/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Systematic Biology Oxford University Press

# Ranked Tree Shapes, Nonrandom Extinctions, and the Loss of Phylogenetic Diversity

, Volume Advance Article – Apr 16, 2018
16 pages

/lp/ou_press/ranked-tree-shapes-nonrandom-extinctions-and-the-loss-of-phylogenetic-IP0z7FnYmt
Publisher
Oxford University Press
ISSN
1063-5157
eISSN
1076-836X
D.O.I.
10.1093/sysbio/syy030
Publisher site
See Article on Publisher Site

### Abstract

Abstract Phylogenetic diversity (PD) is a measure of the evolutionary legacy of a group of species, which can be used to define conservation priorities. It has been shown that an important loss of species diversity can sometimes lead to a much less important loss of PD, depending on the topology of the species tree and on the distribution of its branch lengths. However, the rate of decrease of PD strongly depends on the relative depths of the nodes in the tree and on the order in which species become extinct. We introduce a new, sampling-consistent, three-parameter model generating random trees with covarying topology, clades relative depths, and clades relative extinction risks. This model can be seen as an extension to Aldous’ one parameter splitting model ($$\beta$$, which controls for tree balance) with two additional parameters: a new parameter $$\alpha$$ quantifying the relation between age and richness of subclades, and a parameter $$\eta$$ quantifying the relation between relative abundance and richness of subclades, taken herein as a proxy for overall extinction risk. We show on simulated phylogenies that loss of PD depends on the combined effect of all three parameters, $$\beta$$, $$\alpha,$$ and $$\eta$$. In particular, PD may decrease as fast as species diversity when high extinction risks are clustered within small, old clades, corresponding to a parameter range that we term the “danger zone” ($$\beta<-1$$ or $$\alpha<0$$; $$\eta>1$$). Besides, when high extinction risks are clustered within large clades, the loss of PD can be higher in trees that are more balanced ($$\beta>0$$), in contrast to the predictions of earlier studies based on simpler models. We propose a Monte-Carlo algorithm, tested on simulated data, to infer all three parameters. Applying it to a real data set comprising 120 bird clades (class Aves) with known range sizes, we show that parameter estimates precisely fall close to the danger zone: the combination of their ranking tree shape and nonrandom extinctions risks makes them prone to a sudden collapse of PD. Beta-splitting model, biodiversity, broken stick, field of bullets model, macroevolution, phylogenetic tree, rarefaction, sampling distribution, self-similar fragmentation As it becomes increasingly clear that human activities are causing a major extinction crisis (Leakey and Lewin, 1995; Glavin, 2007; Wake and Vredenburg, 2008; Barnosky et al., 2011), several theoretical studies have aimed at characterizing how the evolutionary legacy of parts of the tree-of-life, and hence also the genetic diversity able to drive future evolution, will decrease in the face of forthcoming extinctions. This evolutionary component of biodiversity can be measured by the phylogenetic diversity (PD), defined as the sum of the branch lengths of the phylogeny spanned by a given set of taxa (Faith, 1992). This metric is increasingly being used to measure biodiversity and to identify conservation strategies (Veron et al., 2015). Nee and May (1997) were the first to formally investigate the expected loss of PD in the face of species extinctions, by simulating species trees using the Kingman coalescent. They found that 80% of the PD can be conserved even when 95% of species are lost. Further studies showed that the loss of PD is in fact much higher when trees are generated through other models of species diversification, such as the Yule or the birth–death models (Morlon et al., 2011b; Mooers et al., 2012; Lambert and Steel, 2013). These models indeed generate longer pendant edges (i.e., branches that lead to the tips), hence lower phylogenetic redundancy, than in the standard Kingman coalescent (used by Nee and May, 1997). However, Nee and May (1997) also showed that PD is very sensitive to the shape of the species tree (also called its “topology”), with extremely unbalanced trees (“caterpilar trees”) losing much more PD than balanced trees (“bush trees”), due to a lack of phylogenetic redundancy (i.e., the presence of recently diverged sister species). Overall, these results highlighted the sensitivity of the loss of PD in response to species extinctions to both edge lengths and tree shape. In this line, we also expect the correlation between the species richness of clades and their relative ages to have a significant impact on the loss of PD (“clade” standing here for any subtree within the full phylogeny). Here the age of a clade, also called “stem age,” denotes the depth (measured from the present) of its root node (i.e., the node where this clade is tied to the rest of the tree). Under random extinction, since smaller clades are more likely to become extinct first, the consequence of their total extinction on PD will depend on the lengths of pendant edges in these clades compared to those in larger clades. Under models with diversification rates that are constant through time and homogeneous across lineages, the time for speciation hypothesis states that the size of clades is correlated to their age (Magallon and Sanderson, 2001), yet several empirical studies found no correlation between the two (Ricklefs 2007; Rabosky et al. 2012; Sánchez-Reyes et al. 2016; but see McPeek and Brown 2007). The effect of such correlation on the loss of PD has not yet been explored, but should be particularly important in unbalanced phylogenetic trees (exhibiting large variation in the species richness of clades), which dominate empirical data (e.g., Guyer and Slowinski, 1991; Heard, 1992; Guyer and Slowinski, 1993; Slowinski and Guyer, 1993; Mooers, 1995; Purvis, 1996; Mooers and Heard, 1997; Blum and François, 2006). Besides, the loss of PD was shown to be influenced by the distribution of extinction risks within species trees. Several studies showed that accounting for realistic scenarios of species extinctions (considering that species with higher extinction risk as per the International Union for Conservation of Nature (IUCN) Red List status are more likely to go extinct first) predicts proportionately higher losses in PD than scenarios with random extinction risks (e.g., and review, Purvis et al., 2000a; von Euler, 2001; Purvis, 2008; Veron et al., 2015). Extinctions may for example be clustered within certain clades (Bennett and Owens, 1997; McKinney, 1997; Russell et al., 1998; Purvis et al., 2000a; Baillie et al., 2004; Bielby et al., 2006; Fritz and Purvis, 2010), correlated to the age of clades (von Euler, 2001; Johnson et al., 2002; Redding and Mooers, 2006), to the species richness of clades (Russell et al., 1998; Hughes, 1999; Purvis et al., 2000a; Schwartz and Simberloff, 2001; von Euler, 2001; Johnson et al., 2002; Lozano and Schwartz, 2005; Vamosi and Wilson, 2008, assuming in some studies a correlation between rarity and extinction risks), or to speciation rates (Heard and Mooers, 2000). In contrast, most of the theoretical analyses of predictions based on model trees (Nee and May, 1997; Mooers et al., 2012; Lambert and Steel, 2013) have all been based so far on the field of bullets model, which considers equal extinction probabilities across species (Raup et al., 1973; Van Valen, 1976; Nee and May, 1997; Vazquez and Gittleman, 1998; but see Heard and Mooers 2000). One can assume extinction events are independent but not identically distributed across species, as considered in the generalized field of bullets model (Faller et al., 2008). In an exchangeable phylogenetic model in which extinction probabilities are themselves random and independent with the same distribution, this would not affect the overall loss of PD (as both models are stochastically equivalent, Lambert and Steel, 2013). However, as stated by Faller et al. (2008), it is essential to explore models that weaken the strong assumption in the (generalized) field of bullets models that extinction events are randomly and independently distributed among the tips of phylogenetic trees. Here, we hence investigate how the loss of PD is influenced by the two abovementioned factors: (i) the ranked shape of the species tree, characterized by the relation between the richness of clades and their age or depth in the tree and (ii) nonrandom extinctions, characterized by the relation between the richness of clades and the extinction risks within them. Here, “ranked shape” refers to the shape of the tree combined with the additional knowledge of relative depths—the order in which nodes appear in the tree, but to the exclusion of the actual divergence times (e.g., Lambert et al., 2017). We introduce a three-parameter model generating random ranked tree shapes endowed with random numbers summing to one at the tips, interpreted as relative abundances (or geographic ranges) of contemporary species. This model can be seen as an extension to Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001) with two additional parameters: a parameter $$\alpha$$ quantifying the relation between the richness of a clade and its relative age (i.e., the rank of appearance of its root node) termed “age-richness index” hereafter, and another parameter $$\eta$$ quantifying the relation between the richness of a clade and its relative abundance or frequency (i.e., the sum of abundances of the species it encompasses divided by the sum of abundances of all extant species in the phylogeny), termed “abundance-richness index” hereafter. When $$\beta=0$$ and $$\alpha=1$$, the ranked shape of the tree is the same as the ranked shape of a standard coalescent tree or of a Yule tree stopped at a fixed time (see Proposition 1 in Supplementary Appendix S1 available on Dryad). We further assume that extinctions of contemporary species occur sequentially in the order of their abundances, starting with the least abundant species, which roughly reduces to the field of bullets model when $$\eta=1$$ (see Proposition 2 in Supplementary Appendix S1 available on Dryad). The parameters of the model are not supposed to map onto biological processes. Our aim is to produce and describe a broad range of ranked tree shapes and extinction risk distributions. The model nevertheless reflects relevant biological patterns. The imbalance and node order of phylogenies may be affected by diversification rates varying with time (Rabosky and Lovette, 2008b; Moen and Morlon, 2014), with species age (Doran et al., 2006; Hagen et al., 2015; Alexander et al., 2015), or among lineages (Cardillo et al., 2005; Maddison et al., 2007; Alfaro et al., 2009; Morlon et al., 2011a). Extinction risk may cluster within the phylogeny if it is correlated with some species characteristic, such as body size (Gaston and Blackburn, 1995; Johnson et al., 2002) or habitat use (Johnson et al., 2002). Finally, tree imbalance, node order and extinction risk clustering are likely to interact if differences in diversification rates and extinction risks are driven by the same trait, or are the result of a common process (e.g., geographic speciation, Pigot et al. 2010). We explore the rate of decrease of PD as species sequentially become extinct, based on simulated data under variation in all three parameters over a significant range of their possible values. Interestingly, the joint variation of the parameter $$\eta$$ with the ranked shape of species trees (set by parameters $$\beta$$ and $$\alpha$$) affects the clustering of extinction risks and the relationship between extinction risks and clade age (determined by the similarity or dissimilarity of the direction of deviations of $$\alpha$$ and $$\eta$$ from 1). Therefore, considering simultaneous variation in $$\beta$$, $$\alpha,$$ and $$\eta$$ allows us to explore the effects on the loss of PD of the different patterns of nonrandom extinctions observed in empirical data. We therefore provide general predictions on the sensitivity of the evolutionary legacy of clades to extinction, as a function of three simple statistics summarizing tree balance, ranked tree shape and the distribution of extinction risks across clades. Following this exploration, we then propose a Monte-Carlo inference algorithm enabling maximum likelihood estimation of the parameters $$\beta$$, $$\alpha,$$ and $$\eta$$ from real data sets. When tested against simulated data, this algorithm performs reasonably well over a wide range of parameter values for phylogenies with 50 tips or more. The estimates of parameters (beta, alpha, eta) on a real data set of bird family phylogenies and their range size distributions finally reveal empirical patterns clustered within a given parameter zone which make these clades particularly prone to strong loss of PD. Materials and Methods Modeling Ranked Tree Shapes The first version of the model we present allows one to generate random ranked tree shapes, that is tree shapes endowed with the additional knowledge of node ranks. Usually, one can generate random ranked tree shapes by time-continuous branching processes stopped at some fixed or random time, where particles are endowed with a heritable trait influencing birth and death rates. In these models, it is generally not possible to characterize the distribution of the tree shape (for an exception, see Sainudiin and Véber, 2016) or to relate it to known distributions whenever it does not have the shape of the Yule tree (i.e., the tree generated by a pure-birth process). Also, since the same trait is usually responsible for both the tree shape and the order of nodes, it is impossible to disentangle the roles of either of these characteristics on the behavior of the tree in the face of current extinctions. Last, these models do not fulfill a very important property called sampling consistency (usually considered in combination with exchangeability, i.e., ecological equivalence between species). This property ensures that one can equivalently draw a random tree with $$n$$ tips from the distribution or draw a tree with $$n+1$$ tips and then remove one tip at random. The model we propose here has two parameters: $$\beta \in (-2,+\infty)$$ (tree balance index) determines the balance of the tree, similarly as in Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001), and $$\alpha \in (-\infty,+\infty)$$ (age-richness index) sets the relation between the species richness of a clade and its relative age (Fig. 2). Figure 1. View largeDownload slide Illustration of the model generating ranked tree shapes. Construction of the ranked shape of a tree containing $$N=5$$ species. (1) Five random marks $$(U_i)_{i \in \lbrace 1,...,5\rbrace}$$ are drawn uniformly in the interval $$[0,1]$$ (the marks on the bottom line). (2) At each time step (time flowing downwards), we randomly select one interval $$X$$, with each interval $$X_j$$ having a weight $$|X_j|^{\alpha}$$ (in black). Then, we draw a random variable $$R$$ in a beta distribution with parameters ($$\beta+1$$, $$\beta+1$$), and split the selected interval $$X$$ into two subintervals, $$X_{\rm{left}}$$ of size $$R|X|$$ and $$X_{\rm{right}}$$ of size $$(1-R)|X|$$ (dark red mark). (3) Repeating this process over time until all intervals $$X_j$$ contain only one mark leads a tree with a ranked shape. Dotted branches correspond to unsampled subtrees (i.e., there is no mark in the corresponding interval). Figure 1. View largeDownload slide Illustration of the model generating ranked tree shapes. Construction of the ranked shape of a tree containing $$N=5$$ species. (1) Five random marks $$(U_i)_{i \in \lbrace 1,...,5\rbrace}$$ are drawn uniformly in the interval $$[0,1]$$ (the marks on the bottom line). (2) At each time step (time flowing downwards), we randomly select one interval $$X$$, with each interval $$X_j$$ having a weight $$|X_j|^{\alpha}$$ (in black). Then, we draw a random variable $$R$$ in a beta distribution with parameters ($$\beta+1$$, $$\beta+1$$), and split the selected interval $$X$$ into two subintervals, $$X_{\rm{left}}$$ of size $$R|X|$$ and $$X_{\rm{right}}$$ of size $$(1-R)|X|$$ (dark red mark). (3) Repeating this process over time until all intervals $$X_j$$ contain only one mark leads a tree with a ranked shape. Dotted branches correspond to unsampled subtrees (i.e., there is no mark in the corresponding interval). Figure 2. View largeDownload slide Phylogenetic trees simulated for different values of $$\beta$$ (tree balance) and $$\alpha$$ (age-richness index). Node depths are set as in a Yule pure-birth process. Parameter values: $$\beta=-1.5$$ (bottom) or 10 (top), $$\alpha=-10$$ (left) or 10 (right), number of species $$N=30$$, approximation parameter (see Supplementary Appendix S1 available on Dryad) $$\varepsilon=0.001$$. Figure 2. View largeDownload slide Phylogenetic trees simulated for different values of $$\beta$$ (tree balance) and $$\alpha$$ (age-richness index). Node depths are set as in a Yule pure-birth process. Parameter values: $$\beta=-1.5$$ (bottom) or 10 (top), $$\alpha=-10$$ (left) or 10 (right), number of species $$N=30$$, approximation parameter (see Supplementary Appendix S1 available on Dryad) $$\varepsilon=0.001$$. The construction of a tree according to this model is done by following the steps indicated hereunder (illustrated on Fig. 1). We start with $$n$$ uniform, independent random variables $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$ in the interval $$[0,1]$$. For each $$i$$, the mark $$U_i$$ is associated to the tip species labeled $$i$$ in the phylogeny. The procedure consists in sequentially partitioning $$[0,1]$$ into a finite subdivision thanks to random variables independent of the marks $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$, until all marks are in distinct components of the partition. At each step, the new point added to the subdivision corresponds to a split event in the tree. In the beginning, there is only one component in the partition (the interval $$[0,1]$$ itself). Each interval $$X$$ of the partition containing at least two marks among the $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$ is given a weight equal to $$|X|^\alpha$$, where $$|X|$$ denotes the width of the interval $$X$$. One of these intervals is selected with a probability proportional to its weight. Draw a random variable $$R$$ in a Beta distribution with parameters $$(\beta+1, \beta+1)$$. The selected interval $$X$$ of width $$|X|$$ is then split into two disjoint subintervals, $$X_{\rm{left}}$$ and $$X_{\rm{right}}$$, with widths $$|X_{\rm{left}}|=R|X|$$ and $$|X_{\rm{right}}|=(1-R)|X|$$. Each subinterval contains a distinct subset of the marks. The marks in the subinterval $$X_{\rm{left}}$$ determine the tips in the left subtree of the phylogeny, and the marks in the subinterval $$X_{\rm{right}}$$ determine the tips in the right subtree. This step is performed even if one subinterval contains no mark among the $$(U_i)_{i \in \lbrace1,\ldots,n\rbrace}$$, which corresponds to a subtree with no sampled species (i.e., species which are seen in the phylogeny). The order in which the splitting subintervals are selected sets the order of branching events (i.e., nodes) in the tree. If no interval contains more than one mark, the process is stopped. Otherwise, go to Step 1. We can relate the tree shape in this model to well-known distributions. Because $$\alpha$$ has no impact on the way we refine the subdivision, the tree shape generated with our model coincides exactly with the tree shape obtained in Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001) with the same parameter $$\beta$$. For small values of $$\beta$$, the intervals are often split close to one of their extremities, and the resulting tree is unbalanced, converging to the perfectly unbalanced “caterpilar” tree as $$\beta\to -2$$. On the contrary, for large values of $$\beta$$, the intervals are often split close to their middle, and the resulting tree is balanced. We stress that unlike most models, $$\alpha$$ can be tuned independently of $$\beta$$, allowing node ranks to vary while keeping the same tree shape. For small values of $$\alpha$$ (in particular $$\alpha<0$$), the smallest subintervals have a higher probability of being selected, so smaller clades tend to be older. On the contrary, for large values of $$\alpha$$, the largest subintervals have a higher probability of being selected, so smaller clades tend to be younger. We notice that as $$\beta$$ gets close to $$-2$$ the effect of $$\alpha$$ vanishes, since most of the time there is only one subinterval containing more than one mark, so only one subinterval to split. In maximally unbalanced tree shape ($$\beta=-2$$), there is only one ranked tree shape and the order of nodes is fixed, so $$\alpha$$ plays no role. As is well-known, the tree obtained with $$\beta=0$$ has the same shape has the tree generated with the Yule process (Yule, 1925) or the Kingman coalescent (Kingman, 1982) after ignoring node ranks (Nee, 2006; Lambert and Stadler, 2013). When $$\alpha=1$$ in addition to $$\beta=0$$, we show in Supplementary Appendix S1 (Proposition 1) available on Dryad at http://dx.doi.org/10.5061/dryad.mv980, that our model generates the same tree shape with node ranks as Yule trees, which is actually known to be the same as the ranked tree shape of the Kingman coalescent tree. The version of the model we present here only allows simulation of trees with $$\beta>-1$$, as the beta distribution is only defined for positive parameter values. Actually, our model coincides with the ranked tree in a self-similar, binary fragmentation with self-similarity index $$\alpha$$ (which is the very age-richness index $$\alpha$$) and with fragmentation measure $$\int_0^1 \delta_{(x,1-x,0,0,\ldots)}\,x^{\beta+1}(1-x)^{\beta+1}dx$$ (as defined in Bertoin, 2002, 2006), which makes sense as soon as $$\beta>-2$$. In Supplementary Appendix S1 (Proposition 3) available on Dryad, we present an algorithm based on fragmentation processes equivalent to that presented above (using one additional approximation parameter $$\varepsilon$$, the maximal frequency of unsampled clades with insignificant richness, consistently set to 0.001). Albeit less intuitive, this method allows us to simulate trees for all $$\beta>-2$$. Last, it is important to notice that our model is both exchangeable and sampling consistent. It is exchangeable because labels can be swapped without changing the distribution of the tree, since marks all have the same distribution. It is sampling consistent because removing tip labeled $$n+1$$ (or any other tip, by exchangeability) amounts to removing mark $$U_{n+1}$$, which does not modify the ranked tree shape obtained from marks $$(U_i)_{i \in \lbrace1,...,n\rbrace}$$. Incorporating Nonrandom Extinctions In order to parametrize the relation between the richness of a clade and its relative abundance, we now add to the model a new parameter $$\eta\ge 0$$ called “abundance-richness index.” Each time an interval $$X$$ is split into two subintervals, $$X_{\rm{left}}$$ and $$X_{\rm{right}}$$ with widths $$|X_{\rm{left}}|=R|X|$$ and $$|X_{\rm{right}}|=(1-R)|X|$$, each of the two subtrees is granted a part of the relative abundance $$A_X$$ of the parental clade equal to   \begin{align*} A_{X_{\rm{left}}} & = \frac{|X_{\rm{left}}|^\eta}{ |X_{\rm{left}}|^\eta+|X_{\rm{right}}|^\eta} \,A_X = \frac{R^\eta}{R^\eta+(1-R)^\eta}\, A_X \\ A_{X_{\rm{right}}} & = \frac{|X_{\rm{right}}|^\eta}{ |X_{\rm{left}}|^\eta+|X_{\rm{right}}|^\eta}\, A_X = \frac{(1-R)^\eta}{R^\eta+(1-R)^\eta} \, A_X. \end{align*} This way of allocating frequencies to taxa is reminiscent of the “broken stick model” (MacArthur, 1957; MacArthur and Wilson, 1967; Colwell and Lees, 2000), where the unit interval is broken into subintervals each representing the frequency or resource share of each species or clade in the community, even though in our case the allocation of relative abundances has no mechanistic interpretation. This is usually done by throwing uniform points independently in the interval or by throwing the points sequentially, always to the right of the last one, leading to the Poisson–Dirichlet distribution appearing in mathematical population genetics (Feng, 2010; Ewens, 2012) as well as in the neutral theory of biodiversity (Hubbell, 2001). The model remains sampling-consistent insofar as each $$A_X$$ is interpreted as the relative abundance of an entire clade, that is the sum of relative abundances of all species, sampled or not (seen or not in the phylogeny), belonging to this clade. Sampling consistency now means that generating a ranked tree shape with relative abundances on $$n$$ tips is equivalent to the following process: generate a ranked tree shape with relative abundances on $$n+1$$ tips, remove one tip at random and sum the abundance of the removed tip to that of its sister clade (i.e., the clade descending from the interior node connected to the removed tip by a pendant edge). In the extinction numerical experiment, we determine the order of species extinctions deterministically based on their abundance: rarest species become extinct first, whereas most abundant species become extinct last. The distribution of relative abundances across species is captured in the abundance-richness index $$\eta$$ (Fig. 3 and in the Supplementary Appendix S1 available on Dryad): when $$\eta=1$$, all tip species have the same relative abundance in expectation ($$\frac{1}{n}$$), such that the correlation between the relative abundance of a clade and its richness would also be close to $$1.0$$. When we simulate extinction with $$\eta=1$$, all species have an approximately equal chance of being removed, and we approach a field of bullets model (see Proposition 2 in Supplementary Appendix S1 available on Dryad; in the case $$\beta=-1$$, the equivalence is exact). For $$\eta>1$$, species in larger clades have larger abundances on average (and thus lower extinction risks), and for $$\eta<1$$, species in larger clades have smaller abundances (and thus higher extinction risks) on average. This modeling approach allows us to tune the sign and strength of the relation between the richness of a clade and the extinction risk of its tip species. Figure 3. View largeDownload slide Distribution of species frequencies across the tips of phylogenetic trees for different values of abundance-richness index $$\eta$$. Dot sizes sort species according to their frequency (larger dots for more abundant species). Parameter values: $$\eta=0.2$$, 1, or 3 (from left to right), $$\beta=0$$, $$\alpha=0$$, number of species $$N=30$$, approximation parameter $$\varepsilon=0.001$$. Results with $$\beta=-1.9$$ are shown in the Supplementary Appendix S1 available on Dryad. Figure 3. View largeDownload slide Distribution of species frequencies across the tips of phylogenetic trees for different values of abundance-richness index $$\eta$$. Dot sizes sort species according to their frequency (larger dots for more abundant species). Parameter values: $$\eta=0.2$$, 1, or 3 (from left to right), $$\beta=0$$, $$\alpha=0$$, number of species $$N=30$$, approximation parameter $$\varepsilon=0.001$$. Results with $$\beta=-1.9$$ are shown in the Supplementary Appendix S1 available on Dryad. Testing the Effect of $$\beta$$, $$\alpha$$, and $$\eta$$ on PD Loss The effect of all three model parameters on the relationship between species loss and PD loss is studied in a systematic way by simulation. We considered values of $$\beta$$ in $$(-2,10]$$, values of $$\alpha$$ in $$[-3,3]$$ and $$\eta$$ in $$[0.1,3]$$. Because our model specifies how interior nodes are ranked in time but not their actual timing, we use a pure-birth process to generate node depths, adding the latter on top of ranked tree shapes. The use of another model for generating node depths leads to qualitatively similar results, albeit quantitatively different (as an illustration, we show results with edge lengths set as in the Kingman coalescent in Supplementary Appendices S4 and S6 available on Dryad). For each set of parameter values, we generated 100 trees with 100 tips ($$N=100$$). We sequentially removed extinct species from these trees (in the order of increasing species abundances, as explained earlier), and computed the remaining PD (sum of all branch lengths; Faith, 1992) for increasing fractions of extinct species. Parameter Inference We inferred the parameters $$\beta$$, $$\alpha,$$ and/or $$\eta$$ from simulated or empirical data sets by maximum likelihood. As is already well-known (Aldous, 1996; Blum and François, 2006; Lambert et al., 2017), the likelihood of a labeled tree shape under Aldous’ $$\beta$$-splitting model is explicit. Since the likelihood of the tree shape under our model is the same as in Aldous’ model (and in particular independent of $$\alpha$$ and $$\eta$$), we can use it to estimate $$\beta$$. In contrast, computing the likelihood of the ranked tree shape requires to follow through time the lengths of all intervals of the partition containing marks, which may decrease without separating marks (unsampled species). Given that the likelihood of the ranked tree (with or without tip abundances) with the additional knowledge of interval lengths is explicit, we use a Monte-Carlo data augmentation procedure, in which the augmentation data are the numbers and sizes of unsampled splits on each branch (which allow us to reconstruct the interval lengths through time). The likelihood of the ranked tree with tip abundances is then computed by averaging over augmentations and is optimized over possible values of $$(\alpha, \eta)$$. We first tested our ability to infer the model parameters on simulated trees. To do so we simulated trees with 20, 50 and 100 tips for all possible combinations of $$\alpha$$ in $$\lbrace -1 , 0, 1, 2 \rbrace$$, $$\beta$$ in $$\lbrace -1 , 0, 1 \rbrace$$ and $$\eta$$ in $$\lbrace 0.2 , 0.5, 1, 1.5 , 2 \rbrace$$. For each tree size and parameter combination, we simulated 20 trees with tip abundances, for a total number of 3600 trees. We then inferred the model parameters on these trees and compared them to the values used in the simulations. The inference of the parameter $$\beta$$ was straightforward, being computed as the maximum likelihood estimate on the interval $$\left[ -2,10\right]$$ with the function maxlik.betasplit from the R-package apTreeshape (Bortolussi et al., 2006). The parameters $$\alpha$$ and $$\eta$$ were estimated with the method introduced hereabove, with values respectively constrained in the intervals $$\left[-4,4 \right]$$ and $$\left[0.1,10 \right]$$. The value of $$\varepsilon$$ (approximation parameter, see Supplementary Appendix S1 available on Dryad) was here again fixed to $$0.001$$. After validating the estimation procedure, we applied it to real bird family trees. We used the Maximum Clade Credibility (MCC) tree from Jetz et al. (2012), and pruned it to keep family level phylogenies. We kept only the phylogenies that included at least 50 species, and used range sizes from Map of Life (https://mol.org/) as proxies for relative abundances. The value of $$\varepsilon$$ and the constraints on parameter ranges were here the same as in the test on simulated phylogenies. The model was coded—and the analyses of phylogenetic trees were performed—using R (R Development Core Team, 2012) and the R packages cubature (Johnson and Narasimhan, 2013), ape (Paradis et al., 2004), sads (Prado et al., 2015), apTreeshape (Bortolussi et al., 2006), and picante (Kembel et al., 2014). The code is available in the R-package apTreeshape (Bortolussi et al., 2006). The list of available functions is given in Supplementary Appendix S10 available on Dryad. Results Influence of Ranked Tree Shape on PD Loss Here we only address the influence of $$\alpha$$ on PD loss, assuming a field of bullets model for species extinctions ($$\eta=1$$). The expected PD loss is then a convex function of the fraction $$p$$ of extinct species (as proved mathematically for any binary tree under the field of bullets model, see Eq (34) in Lambert and Steel, 2013), always lying below $$p$$ (Fig. 4A,C,E,G). Figure 4. View largeDownload slide Influence of the ranked tree shape (tree balance $$\beta$$ and age-richness index $$\alpha$$) on PD loss, for increasing fractions of species extinctions. Tree balance $$\beta$$ changes from $$10$$ (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”). Results are shown either as a function of the extinction fraction $$p$$ (left column; for different $$\alpha$$ values) or as a function of $$\alpha$$ (right column; for different extinction fractions $$p$$). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right in a, c, e, g; from light to dark color in b, d, f, h). The dotted lines in a, c, e, g show the bisector. Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 4. View largeDownload slide Influence of the ranked tree shape (tree balance $$\beta$$ and age-richness index $$\alpha$$) on PD loss, for increasing fractions of species extinctions. Tree balance $$\beta$$ changes from $$10$$ (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”). Results are shown either as a function of the extinction fraction $$p$$ (left column; for different $$\alpha$$ values) or as a function of $$\alpha$$ (right column; for different extinction fractions $$p$$). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right in a, c, e, g; from light to dark color in b, d, f, h). The dotted lines in a, c, e, g show the bisector. Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Consistently with previous studies (Nee and May, 1997; von Euler, 2001) we find that when the relation between age and richness of clades is similar as that in Yule trees ($$\alpha=1$$), very unbalanced trees (caterpilar-like trees) lose more PD in the face of species extinctions than Yule or more balanced trees (Fig. 4G,H vs. A–D, with $$\alpha=1$$). The effect is nonlinear in $$\beta$$: the tree shape has little influence on the loss of PD when $$\beta\ge -1$$, but increases sharply as $$\beta$$ decreases from $$-1$$ to $$-1.9$$ (results as a function of $$\beta$$ in Supplementary Appendix S2 available on Dryad). Unbalanced tree shapes are associated with the presence of long edges leading to evolutionary distinct species (Fig. 2). These edges constitute an important fraction of the PD in unbalanced species trees, so that their extinction generates a significant drop in PD. As $$\beta$$ gets closer to $$-2$$ (case of the “caterpilar tree”), the expected PD loss approaches the fraction of extinct species (Fig. 4G). Considering ranked tree shapes shows, however, that the order of nodes has a significant influence on the loss of PD, and on the effect of $$\beta$$ on this loss. If the age and richness of clades are positively correlated ($$\alpha>0$$), the loss of PD is reduced, especially at intermediate extinction fractions (Fig. 4A–F). This is because the smallest subtrees, more prone to early extinction, are younger and hence contain a lower fraction of the PD (Fig. 2). If the age and richness of clades are negatively correlated ($$\alpha<0$$), the loss of PD rises, especially at intermediate extinction fractions. The smallest subtrees, prone to extinction, are older and hence contain more evolutionary distinct species (Fig. 2). This generates losses of PD similar to those observed when the tree shapes are very unbalanced (PD loss equal to the fraction of extinct species). As expected, the effect of $$\alpha$$ is evened out in very unbalanced trees ($$\beta$$ close to $$-2$$; Fig. 4G,H), for which the loss of PD remains close to its highest value whatever the value of $$\alpha$$. In the case of the maximally unbalanced tree shape, there is only one ranked tree shape and the order of nodes is fixed. All these effects of ranked tree shapes on the loss of PD are qualitatively conserved if node depths are distributed as in the Kingman coalescent (instead of the Yule process). In the case of Yule trees, PD loss slightly increases with the initial size of the tree, an effect which is due to more efficient sampling of large values in the common (exponential) distribution of node depths. Yet the results presented above are qualitatively conserved if the size of phylogenetic trees changes (analyses performed with number of species $$N=50$$ and $$N=200$$; see the Supplementary Appendices 3 and 4 available on Dryad). Influence of Nonrandom Extinction Risks on PD Loss The strength of the relation parameterized by $$\eta$$, between the richness of a clade and its relative abundance (here directly influencing the extinction risk of its species) may have a paramount influence on the loss of PD in the face of extinctions (Fig. 6). In trees with ranked tree shapes similar to Yule trees ($$\beta=0$$, $$\alpha=1$$), the concentration of high extinction risks in small clades ($$\eta>1$$) increases the loss of PD, by promoting the extinction of entire clades (Fig. 5). In contrast, when extinction risks are higher in species-richer clades ($$\eta<1$$), phylogenetic redundancy (and hence the likelihood of conserving at least one species per subtree) limits the loss of PD until high extinction levels. Figure 5. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss in Yule trees, for increasing fractions of species extinctions $$p$$. Results are shown either a) as a function of the extinction fraction $$p$$ (for different $$\eta$$ values, with dotted lines showing the bisector) or b) as a function of $$\eta$$ (for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Parameter values: $$\beta=0$$, $$\alpha=0$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 5. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss in Yule trees, for increasing fractions of species extinctions $$p$$. Results are shown either a) as a function of the extinction fraction $$p$$ (for different $$\eta$$ values, with dotted lines showing the bisector) or b) as a function of $$\eta$$ (for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Parameter values: $$\beta=0$$, $$\alpha=0$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 6. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss, for different ranked tree shapes and increasing fractions of species extinctions. Tree balance $$\beta$$ ranges from 10 (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”), and age-richness index $$\alpha$$ ranges from $$-2$$ (A, C, E, G) to $$2$$ (B, D, F, H). Results are shown either as a function of the extinction fraction $$p$$ (left side; for different $$\eta$$ values, and with dotted lines showing the bisector) or as a function of $$\eta$$ (right side; for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 6. View largeDownload slide Effect of abundance-richness index $$\eta$$ on PD loss, for different ranked tree shapes and increasing fractions of species extinctions. Tree balance $$\beta$$ ranges from 10 (top row, “bush trees”) to $$-1.9$$ (bottom row, “caterpilar trees”), and age-richness index $$\alpha$$ ranges from $$-2$$ (A, C, E, G) to $$2$$ (B, D, F, H). Results are shown either as a function of the extinction fraction $$p$$ (left side; for different $$\eta$$ values, and with dotted lines showing the bisector) or as a function of $$\eta$$ (right side; for extinction fractions $$p$$ increasing from 0.01 to 0.98 from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. The effect of $$\eta$$ is modified by the ranked shape of species trees. The strength of the age-richness relation (set by $$\alpha$$) modulates the additional loss of PD induced by $$\eta>1$$ (i.e., lower abundances in smaller clades; Fig. 6A–F). When $$\alpha<0$$, smaller clades are not only more prone to extinction but also have deeper nodes, hence more evolutionary distinct species, which increases even further the loss of PD. Unlike in the field of bullets model, the expected PD loss as a function of the fraction $$p$$ of extinct species can even change from convex to concave, and so take values larger than $$p$$ (Fig. 6C,E). When $$\alpha>0$$, smaller clades are more prone to extinction but have shallower nodes, which counteracts the increase of PD loss due to $$\eta>1$$. To summarize, PD loss is increased when $$\eta>1$$ compared to $$\eta=1$$, with a maximal effect for negative values of $$\alpha$$, progressively flattening as $$\alpha$$ grows. We call the “danger zone” the region of parameters corresponding to the theoretical phylogenies that suffer a maximal rate of PD loss straight from the first few extinction events, that is, close to 1% of PD lost for the first 1% of species lost. In the plane $$(\alpha, \eta)$$, the “danger zone” corresponds to $$\{\alpha<0, \eta >1\}$$. As testified by Fig. 6, phylogenies in this zone can even suffer a rate of PD loss which is larger than 1:1 from the first extinction and sustains itself above 1:1 throughout the extinction crisis. In contrast, $$\alpha$$ has little effect on the decrease in PD loss induced by $$\eta<1$$ (i.e., higher abundances in small clades). Indeed, when $$\eta<1$$, the deepest nodes are always protected regardless of the value of $$\alpha$$: when $$\alpha<0$$ the deepest nodes are in small clades which are protected from extinctions by the high relative abundances of its species (due to $$\eta<1$$); when $$\alpha>0$$, the deepest nodes are in large clades which are protected by phylogenetic redundancy. The influence of $$\eta$$ on PD loss is amplified by unbalanced tree shapes ($$\beta<0$$; Fig. 6E–H) and buffered by balanced tree shapes ($$\beta>0$$; Fig. 6A,B), because lower values of $$\beta$$ enhance richness inequalities between clades and raise in turn the influence of $$\eta$$ on PD loss. This interaction between parameters $$\eta$$ and $$\beta$$ overwhelms the influence of $$\alpha$$ (Fig. 6). In the plane $$(\beta, \eta)$$, the “danger zone” is $$\{\beta<-1, \eta >1\}$$ and the previous remark thus implies that in the 3D parameter space, the danger zone is $$\{\alpha <0$$ or $$\beta<-1; \eta >1\}$$. Interestingly, the effect of $$\beta$$ is highly dependent on how extinction risks are distributed within the phylogeny (Fig. 7, and results with other $$\alpha$$ values in Supplementary Appendix S7 available on Dryad). For $$\eta=1$$, we recover the well-known pattern of decreased PD loss as the tree gets more balanced. However, for $$\eta<1$$ we see the reverse pattern, that is PD loss increases with the balance of the tree. Recall that $$\eta <1$$ buffers PD loss, because extinction risks are clustered in the species-richer clades which also display higher phylogenetic redundancy (smaller pendant edges). When the tree is maximally unbalanced, $$\eta<1$$ causes the longest pendant edge to subtend the tip with the largest abundance (and hence to be the last to become extinct). Therefore, the order of extinctions coincides exactly with the increasing order of pendant edge lengths, which results in minimal PD loss for any given level of extinction. In a more balanced phylogeny, the distribution of clade sizes is more even and the buffering effect of the clustered extinction on PD loss is reduced. Figure 7. View largeDownload slide Effect of tree balance $$\beta$$ on PD loss, for different values of clades abundance-richness indexes $$\eta$$, and increasing fractions of species extinctions $$p$$. The clades abundance-richness index $$\eta$$ ranges from $$0.2$$ (left) to 3 (right), and the extinction fraction $$p$$ increases from 0.01 to 0.98 (from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: clades age-richness index $$\alpha=-2$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. Figure 7. View largeDownload slide Effect of tree balance $$\beta$$ on PD loss, for different values of clades abundance-richness indexes $$\eta$$, and increasing fractions of species extinctions $$p$$. The clades abundance-richness index $$\eta$$ ranges from $$0.2$$ (left) to 3 (right), and the extinction fraction $$p$$ increases from 0.01 to 0.98 (from light to dark color). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: clades age-richness index $$\alpha=-2$$, number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$. For $$\eta >1,$$ we again recover the well-known pattern of decreased PD loss with increasing $$\beta$$. However, when we also have $$\alpha<0$$, the relationship between PD loss and $$\beta$$ is not monotonic, that is for any particular level of extinction, the maximal PD loss is reached for trees with intermediate balance. Recall that $$\alpha<0$$ causes small clades to be relatively older and so to contribute more to PD. The maximal loss of PD thus occurs when extinction risks cluster in small clades. And indeed, when $$\eta>1$$, at each splitting event the species-richer subtree gets a bigger abundance than the species-poorer subtree. However, within a given clade, the abundance of a species should decrease with the number of nodes (splitting events) on its lineage. This latter effect is stronger in unbalanced trees; in balanced trees, extinction risks cannot cluster in small clades, due to the absence of small clades. Trees with intermediate balance do display small clades, and these small clades are large enough to share their low abundance ($$\eta >1$$) into a few species with very low abundance. These species go extinct first, resulting in maximal PD loss. Effect of Species Extinctions on Tree Shape We study the effect of species extinctions on tree shape, seeking in particular to check if the influence of $$\eta$$ on the patterns of PD loss can be explained by changes in tree shape as species become extinct. Figure 8 shows the balance (defined here as the maximum likelihood estimate $$\hat\beta$$ of the parameter $$\beta$$) of the species tree estimated after a fraction $$p$$ of its species have become extinct. When $$\eta=1$$, tree balance is very little altered by extinctions except in very balanced trees, as predicted by the sampling consistency of the model ($$\eta=1$$ amounts to removing species at random except when $$\beta\gg 1$$, see Supplementary Appendix S1 available on Dryad). When $$\eta<1$$, trees tend to become more and more balanced as $$p$$ increases ($$\hat\beta$$ increases with $$p$$), whereas when $$\eta>1$$ trees tend to become more and more similar to Yule trees ($$\hat\beta\to 0$$ as $$p\to 1$$). The effect of $$\eta$$ on PD loss cannot be reduced to its effect on changes in tree shape due to extinctions. On the one hand, $$\eta$$ mostly affects the shape of trees with $$\beta>-1$$ (Fig. 8), whereas tree shape has most effect on PD loss when $$\beta$$ varies between $$-2$$ and $$-1$$ (Fig. 4A,C,E with $$\alpha=0$$). In addition, if the effect of $$\eta$$ on tree shape had a significant influence on PD loss, $$\eta>$$1 should increase this loss when $$\beta>0$$ (by decreasing the balance of trees; Fig. 8D) and decrease it when $$\beta<0$$ (by increasing the balance of trees). Yet, the changes we observe in the effect of $$\eta>1$$ on PD loss for different $$\beta$$ values are the reverse of this prediction. Therefore, the indirect effects of $$\eta$$ (through changes in tree shape) are negligible compared to its direct effects (through nonrandom distribution of extinction risks). Figure 8. View largeDownload slide Effect of abundance-richness index $$\eta$$ on the balance of phylogenetic trees after extinctions (Maximum Likelihood Estimate $$\hat\beta$$ of $$\beta$$). Initial tree balance $$\beta$$ ranges from 10 (brown dots and lines, “bush trees”) to $$-1.9$$ (green dots and lines, “caterpilar trees”). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$, $$\alpha=0$$. Figure 8. View largeDownload slide Effect of abundance-richness index $$\eta$$ on the balance of phylogenetic trees after extinctions (Maximum Likelihood Estimate $$\hat\beta$$ of $$\beta$$). Initial tree balance $$\beta$$ ranges from 10 (brown dots and lines, “bush trees”) to $$-1.9$$ (green dots and lines, “caterpilar trees”). Extinction fraction $$p$$ increases from 0.01 to 0.98 (from left to right). Results are based on 100 simulation replicates: plain lines give median values and light areas give 95% confidence intervals. Other parameter values: number of species $$N=100$$, approximation parameter $$\varepsilon=0.001$$, $$\alpha=0$$. As reported above results on the effects of nonrandom extinctions on the loss of PD are conserved when node depths are distributed as in the Kingman coalescent, or when the size of phylogenetic trees changes (analyses performed with $$N=50$$ and $$N=200$$; see Supplementary Appendices S5 and S6 available on Dryad). Parameter Inference When tested against simulated data, the Monte-Carlo inference algorithm by data augmentation performs reasonably well on phylogenies with more than 50 tips for a wide range of parameters (see Supplementary Appendix S8 available on Dryad). As expected, the estimation of $$\beta$$ on trees with at least 50 tips is accurate, since the likelihood formula of the unranked tree is explicit, and this accuracy increases as $$\beta$$ decreases. The inference algorithm also returns overall good estimates of $$\eta$$ and $$\alpha$$ whenever $$\eta>0.3$$. The inference of $$\alpha$$ is unbiased except in the cases where $$\beta < 0$$ and $$\eta < 0.3$$. This corresponds to cases where the unsampled clades are numerous because $$\beta$$ is small, and they have a strong impact on the reconstruction of intervals because $$\eta$$ is small. The inferred $$\eta$$ is overestimated for trees with only $$50$$ tips. For $$\beta < 0$$ and $$\alpha \ge 0$$, $$\eta$$ is slightly overestimated whatever the tip number. For $$\beta > 0$$ and $$\alpha \leq 0,$$ inferences are good for trees with at least $$100$$ tips. Empirical Values Estimates of parameter values on real data shows consistent patterns across all bird family trees. Unsurprisingly, we find negative $$\beta$$ values, mostly comprised between 0 and $$-1$$, corresponding to unbalanced trees (see Supplementary Appendix S9 available on Dryad). Since the estimation of $$\beta$$ is quite accurate for low true values of $$\beta$$ and is biased towards larger estimates than the true value otherwise, these estimates can be taken with confidence. The estimates of $$\eta$$ vary between $$1$$ and $$1.5$$. This indicates that, within bird families, species in small clades tend to have smaller range sizes than species in larger clades. The above study showed that low $$\eta$$ values can be difficult to detect in unbalanced trees. Yet when this is the case, $$\eta$$ is found to be close to the maximal value allowed in the inference (here $$10$$), which is not the case here. We can therefore be confident that these values do not reflect a bias in the inference, but reflect a true pattern in the distribution of range sizes within the phylogenies. Finally, the estimates of $$\alpha$$ are clustered around $$0$$, indicating that there is no correlation between clade sizes and clade depths within each bird phylogeny. This in contrast with what is expected in most explicit models of diversification, where larger clades take more time to diversify, resulting in a strong positive correlation between the depth and the size of clades. When jointly inferring $$\alpha$$ and $$\eta$$ the choice to use range size to infer $$\eta$$ is likely to have an impact on the inferred $$\alpha$$ (because the values of the intervals are reconstructed using tip values, inappropriate tip values would lead to incorrect $$\alpha$$). Therefore, we also ran the inference of $$\alpha$$: we find fairly similar results between values obtained with the inference of $$\alpha$$ only compared to the full inference (the median of the inferred $$\alpha$$ for trees with at least $$50$$ tips is $$0.19$$ when $$\alpha$$ is inferred alone and $$0.05$$ when both $$\alpha$$ and $$\eta$$ are inferred), indicating that tree shape is indeed driving the result (Supplementary Appendix S9 available on Dryad equivalent as Fig. 9 with the $$\alpha$$ inferred without knowledge of the tip range sizes). Figure 9. View largeDownload slide Inferred model parameters on bird family trees of 50 tips or more. Maximum posterior estimate of age-richness index $$\alpha$$ (x-axis), maximum posterior estimate of abundance-richness index $$\eta$$ (y-axis) and maximum likelihood estimate of tree balance $$\beta$$ (point color). Point sizes are proportional to the number of tips in trees, $$N$$. The dashed vertical line shows the value of $$\alpha$$ for trees generated by a birth–death model, and the dashed horizontal line shows the value of $$\eta$$ for which extinction probabilities are distributed within the tree as in a field of bullets model. For all inferences, approximation parameter $$\varepsilon$$ was set to $$0.001$$. Figure 9. View largeDownload slide Inferred model parameters on bird family trees of 50 tips or more. Maximum posterior estimate of age-richness index $$\alpha$$ (x-axis), maximum posterior estimate of abundance-richness index $$\eta$$ (y-axis) and maximum likelihood estimate of tree balance $$\beta$$ (point color). Point sizes are proportional to the number of tips in trees, $$N$$. The dashed vertical line shows the value of $$\alpha$$ for trees generated by a birth–death model, and the dashed horizontal line shows the value of $$\eta$$ for which extinction probabilities are distributed within the tree as in a field of bullets model. For all inferences, approximation parameter $$\varepsilon$$ was set to $$0.001$$. Discussion A New Integrative Measure of the Age-Richness Relation, $$\alpha$$ We introduced here a new model for random ranked tree shapes with a fixed, arbitrary number of tips. This model features two parameters, $$\beta$$ and $$\alpha$$ tuning respectively the balance of the tree and the relative depths or ages of its nodes. Trees with $$\beta\le 0$$ are unbalanced and trees with $$\beta>0$$ are balanced. Whatever the value of $$\alpha$$, the shape of the tree is the same as in Aldous’ $$\beta$$-splitting model (Aldous, 1996, 2001). Large clades coalesce deep in the tree when $$\alpha>0$$ and are shallower than smaller clades when $$\alpha<0$$. When $$\beta=0$$ and $$\alpha=1$$, the tree has the same ranked shape as the Kingman coalescent and the Yule tree. Our model does not explicitly incorporate biological processes but enables to generate a broad range of tree shapes by decoupling tree balance and the relation between age and size of clades. In addition, this model is the first model (except the two aforementioned models and the trivial case of the “caterpilar tree”) for ranked tree shapes satisfying sampling-consistency, in the sense that a tree with $$n$$ tips has the same distribution as a tree with $$n+1$$ tips with one tip removed at random. This property is essential to ensure the robustness of the model with respect to incomplete taxon sampling (Heath et al., 2008; Cusimano et al., 2012; Stadler, 2013). Predictions from this model highlight the importance of accounting for node ranks to understand forthcoming changes in macroevolutionary patterns of PD. They show in particular that the relationship between the species richness of a clade and its relative depth in the tree, set by parameter $$\alpha$$ in the model, can have profound impacts on the rate of PD loss (Fig. 4). This parameter $$\alpha$$ constitutes a new index quantifying the age-richness relation of subclades within a phylogeny. A large number of studies have already considered the age-size correlation, assessing its existence (significance, sign, and pattern) across multiple phylogenetic trees, based on one value of species richness and crown or stem age per phylogeny (e.g., Magallon and Sanderson, 2001; Bokma, 2003; Ricklefs, 2006; McPeek and Brown, 2007; Rabosky et al., 2012; Sánchez-Reyes et al., 2016). These studies notably aimed at testing the hypothesis of time-limited diversity patterns, versus hypotheses of diversity set by diversification rates or by limits to diversity (McPeek and Brown, 2007; Ricklefs, 2007; Rabosky, 2009; Barraclough, 2010; Sánchez-Reyes et al., 2016). Our new index $$\alpha$$ is different in that it can be measured by maximizing the likelihood on a single phylogeny, implicitly integrating over all subclades of this phylogeny. An interesting consequence is that one does not have to choose which clades to include in the analysis. For example, $$\alpha$$ is not sensitive to the definition of higher taxa (Stadler et al., 2014; Sánchez-Reyes et al., 2016). Moreover, similarly to what can be done with the index $$\beta$$ (compared to other measures of tree balance; Kirkpatrick and Slatkin, 1993; Aldous, 1996, 2001), a measure of the age-richness relation in a given phylogeny is provided by the maximum likelihood estimate of the model-based parameter $$\alpha$$. Last, we stress that our model does not require the precise knowledge of node datings in the phylogeny but only the relative positions of nodes in time, which preserves $$\alpha$$ estimates from the inaccuracies of time calibrations (Kumar, 2005; Welch and Bromham, 2005; Pulquério and Nichols, 2007; Forest, 2009; Schwartz and Mueller, 2010). Ranked Tree Shapes and the Loss of PD Our results confirm that in the field of bullets model very unbalanced trees undergo stronger loss of PD than balanced trees, under equal fraction of species extinctions. This property was already well known (Nee and May, 1997) but is important to recall given the predominance of unbalanced phylogenetic trees in nature ($$\beta$$ values being often close to $$-1$$; e.g., Guyer and Slowinski, 1991; Heard, 1992; Guyer and Slowinski, 1993; Slowinski and Guyer, 1993; Mooers, 1995; Purvis, 1996; Mooers and Heard, 1997; Blum and François, 2006). However, our results also show that the temporal order of nodes among subtrees (set by the parameter $$\alpha$$) may have even stronger effects than tree balance (set by parameter $$\beta$$; compare the effect of the latter in Supplementary Appendix S8 available on Dryad to that of $$\alpha$$ on Fig. 4). Besides, $$\alpha$$ values below $$0$$ cause drops of PD almost as abrupt as those observed with “caterpilar” shapes ($$\beta$$ close to $$-2$$, with $$\alpha=1$$; Fig. 4D,H). It is therefore essential to consider the ranked shapes of species trees to understand the expected patterns of loss of PD. Values of $$\alpha$$ deviating from $$1$$ may arise from differences in stages of diversification among subtrees, resulting from heterogeneity in biotic or abiotic factors acting on diversification processes in different parts of the species tree. This could be due to bursts of diversification in certain subtrees (e.g., following from key innovations or from migration to empty spatial or ecological space), either recently (resulting in $$\alpha<0$$) or early in the history of clades (resulting in $$\alpha>0$$). Alternatively, $$\alpha$$ values deviating from $$1$$ could be linked to changes in extinction rates in distinct parts of the tree (e.g., due to changes in the biotic or abiotic environment of phylogenetically related species sharing similar ecological niches). Age-dependent speciation (Venditti et al., 2010; Hagen et al., 2015) and extinction (Pearson, 1992; Alexander et al., 2015) are also likely to make node ranking deviate from what is expected in a homogeneous birth–death model. Heterogeneity in diversification rates across the species tree associated with asymmetric competition among species (e.g., evolutionary advantage to previously established species) could limit diversification in younger subtrees, hence leading to $$\alpha>0$$. Last, $$\alpha$$ can be found negative due to the presence of relictual lineages, that is, old clades harboring few species surviving to the present. Modeling Nonrandom Extinctions: $$\eta$$ and the Loss of PD The incorporation of the abundance-richness index $$\eta$$ within the framework provided by Aldous’ $$\beta$$-splitting model allowed us to go beyond the field of bullets assumption. The model allows one to simulate a trait covarying with the shape of the phylogeny. We have so far interpreted this trait as species relative abundance or relative range size, but it could be any species trait with values summing to $$1$$ at the level of the phylogeny. In particular, the trait value of a species can be interpreted as a probability of sampling this species, which is consistent with the initial interpretation of the trait as a relative abundance. In passing, we devised a model of abundance distributions (equivalently interpreted as range size distributions) covarying with the phylogeny, in the broken-stick tradition (MacArthur, 1957; MacArthur and Wilson, 1967). When $$\eta>1$$, the most abundant species are in species-rich clades whereas when $$\eta<1$$ the most abundant species are in species-poor clades. When $$\eta=1$$ all species have the same abundance on average. Here, extinctions are assumed to occur sequentially in the order of increasing abundances. By doing this, we only consider the part of extinction risks that is due to the rarity of a species. In nature, relative extinction risk indeed depends on species abundance, but also on many other features (e.g., dynamics of population growth or decline, fragmentation into subpopulations, biotic or abiotic changes; IUCN, 2012), and may have a significant stochastic component. The simple framework we use to determine extinctions allows us to focus on the direct impact of the distribution of ranked abundances within trees on the loss of PD. This framework can easily be modified to include extrinsic causes of extinctions by taking extinction risk to be a function of abundance and additional known factors. Previous studies concluded that PD loss is increased if extinction risks are clustered in the phylogeny (Davies and Yessoufou, 2013), but that this effect is not substantial (Parhar and Mooers, 2011). Our model shows that the effect on PD loss depends on the way these extinction risks are distributed among clades: PD loss is increased by $$\eta>1$$ (i.e., higher extinction risks in small clades; Fig. 6). Such a distribution of extinction risks may arise from subtrees having low species richness because of higher extinction rates, either due to intrinsic factors (species features that would make them more susceptible to extinction; e.g., long generation time, or low variance or phenotypic plasticity of key ecological traits providing resistance to perturbations or evolutionary advantages in relation to biotic interactions; Purvis et al., 2000c; Johnson et al., 2002) or to extrinsic factors (threats affecting the spatial or ecological space shared by species of the subtree; e.g., Russell et al., 1998; Hughes, 1999; Purvis et al., 2000c; von Euler, 2001; Johnson et al., 2002). Higher extinction risks in small subtrees could also be due to resource limitation affecting simultaneously the density of individuals and the diversity of species, and hence demographic stochasticity (the island effect, Mooers et al., 2009); or to stabilizing selection (e.g., due to competition or to the absence of available spatial or ecological space in the surrounding environment), limiting adaptation and increasing species vulnerability in the face of perturbations (Purvis et al., 2000c; Purvis, 2008). In contrast, $$\eta<1$$ buffers the loss of PD. Higher extinction risks in larger subtrees could result from a trade-off between species richness and average species abundance, provided constrained metacommunity size (with variation along this trade-off following for instance from landscape structure and dynamics, such as geographical isolation affecting the occurrence of allopatric speciation events), from recent speciation events associated with a decrease in average species abundance, geographical range or niche width, or from recent extinction events that removed the most extinction-prone species from certain clades (leaving the latter smaller and with less extinction-prone species; Schwartz and Simberloff, 2001; Lozano and Schwartz, 2005). Hence, $$\eta$$ is expected to vary across clades according to the metacommunity structure and the underpinning diversification dynamics. Given its striking effects on PD loss, this factor should also be accounted for to understand potential future losses of PD. Combined Effects of $$\beta$$, $$\alpha$$, and $$\eta$$: Reversing Some Expected Patterns of PD Loss The influence of $$\eta$$ on the loss of PD is enhanced by $$\alpha<0$$ (small clades containing evolutionary distinct species) and $$\beta<0$$ (more variability in clades richness) (Fig. 6). However, a stronger clustering of extinction risks does not necessarily lead to higher loss of PD (e.g., if extinctions occur first in species-richer subtrees—which contain more phylogenetic redundancy—as in the case when $$\alpha>0$$ and $$\eta<1$$). These interactions between the effects of $$\beta$$, $$\alpha,$$ and $$\eta$$ may reverse two well-known patterns of variation in the loss of PD (Nee and May, 1997). First, the increase in PD loss with tree imbalance can be hampered by $$\eta$$ values deviating from one (Fig. 7 and Supplementary Appendix S7 available on Dryad). In particular when $$\eta<1$$, this pattern results from the preferential extinction of phylogenetically redundant species in more unbalanced trees when extinction risks are clustered in large clades. Second, when $$\eta>1$$ and $$\alpha<0$$ the loss of PD proceeds faster than that of species diversity (turning their relationship from convex to concave, except in very balanced or very unbalanced trees; Fig. 6C,E). This pattern (also highlighted in Heard and Mooers, 2000) is caused by the preferential extinction in small subtrees containing evolutionary distinct species. The only other cases where such high loss of PD is reached is when $$\beta<-1$$ and $$\eta>1$$. This parameter zone is a macroevolutionary danger zone, where particular phylogenetic shapes combine with clustered patterns of extinction probabilities to produce large and rapid losses of evolutionary history, much like ice and wind bringing down thick branches of trees in winter. In this region of the parameter space ($$\beta <-1$$ or $$\alpha <0$$; $$\eta>1$$), phylogenies are prone to a sudden collapse of PD. Loss of PD in Bird Family Phylogenies Our inference study shows that the phylogeny of bird families tend to exhibit $$\beta$$ values comprised between $$-1$$ and $$0$$. A similar result was found in many macroevolutionary studies, commonly observing values of $$\beta$$ clustering around $$-1$$ in real phylogenies (e.g., Guyer and Slowinski, 1991; Heard, 1992; Guyer and Slowinski, 1993; Slowinski and Guyer, 1993; Mooers, 1995; Purvis, 1996; Mooers and Heard, 1997; Blum and François, 2006). With these topologies, we expect both $$\alpha$$ and $$\eta$$ to play a major role in determining the potential losses of PD (Fig. 6E,F). We observed $$\alpha$$ values clustering around zero, consistently with several empirical studies that found no positive relation between clade age and clade size (Ricklefs, 2007, 2009; Rabosky et al., 2012). These values contrast with the value of $$1$$ expected in Yule trees, and make phylogenies very sensitive to PD loss. Our estimates of $$\eta$$ values, based on the distribution of range sizes in bird family phylogenies, all fall between $$1$$ and $$1.5$$. This indicates that species in small clades tend to have smaller ranges than species in bigger clades. Range size has been shown to be one of the most important correlates of extinction risks and is one of the IUCN red list criteria (Purvis et al., 2000b; Cardillo et al., 2006; Lee and Jetz, 2011; IUCN, 2012; Arbetman et al., 2017). Expectedly, range sizes in our data set are significantly different among IUCN status classes (see Supplementary Appendix S9 available on Dryad). This pattern of extinction risks clustering in species-poor clades, highlighted in our data by $$\eta$$ values above $$1$$, has also been shown in plants (Vamosi and Wilson, 2008), and in birds and mammals for past extinctions (Russell et al., 1998; Mooers et al., 2009). Considering the three parameters together, we find that bird family trees are situated close to the region of the parameter space termed “danger zone,” for which we find the loss of PD to be at least as fast as the loss of species diversity. In particular, the combination of negative $$\alpha$$ values with $$\eta>1$$ leads to higher extinction risks for evolutionary distinct species. We can expect such a pattern as a result from evolutionary mechanisms acting simultaneously on different features of trees. For example, subtree-specific susceptibility to extinction, or stabilizing selection generating relictual lineages, are both expected to beget small subtrees with high divergence times also endowed with high species extinction risks. This pattern has been already found for past extinctions in birds using a species level measure of evolutionary distinctiveness; authors observed in that case a similar loss of species and PD (von Euler, 2001; Szabo et al., 2012). Evolutionary distinct bird lineages were also shown to be more threatened by agricultural expansion and intensification than more recent lineages in Costa Rica (Frishkoff et al., 2014). This was also found in other taxa, such as marsupial mammals (Johnson et al., 2002) and Sebastes (Magnuson-Ford et al., 2009) (but see Davies et al. 2011, who show that in South African plants, extinction risks cluster in younger, fast-diversifying genera). A striking result of our inference study relates to the narrow range of $$\alpha$$ values obtained as soon as the trees are large enough for the inference to be accurate (see Supplementary Appendix S9 available on Dryad for inferred parameter values as a function of the tip number in the phylogenies). This value, which differs from what is found in birth–death models, adds a new conundrum concerning the shape of empirical trees. Branch Lengths in Empirical Phylogenies The parameter $$\alpha$$ of the model shapes the order in which speciations take place, but does not instantiate the actual times between two consecutive speciation events, that is, edge lengths. In the numerical investigations of PD loss, we considered two models for edge lengths: the pure-birth process (Yule, 1925) and the Kingman coalescent (Kingman, 1982). Using either of these models did not affect our results qualitatively, but affected them quantitatively (compare Figs. 4 and 6 to figures provided in Supplementary Appendices S4 and S6 available on Dryad). Our modeling framework allows easy exploration of predictions under different models of edge lengths. This is interesting as many empirical phylogenies are not time-calibrated, or imprecisely. Besides, empirical phylogenetic trees were shown to often exhibit a decrease in branching tempo, that is, in the rate of lineage accumulation through time (characterized in particular by estimates of the statistic $$\gamma<0$$; e.g., Nee et al., 1992; Zink and Slowinski, 1995; Lovette and Bermingham, 1999; Pybus and Harvey, 2000; Rüber and Zardoya, 2005; Kozak et al., 2006; Seehausen, 2006; Weir, 2006; McPeek, 2008; Phillimore and Price, 2008; Rabosky and Lovette, 2008a; Jønsson et al., 2012; Moen and Morlon, 2014). Hence, quantitative predictions on the loss of phylogenetic diversity in the face of species extinctions could be further increased by accounting for real branch lengths. Moreover, several theoretical studies suggested that the branching tempo of species trees may change with clade age, decreasing in particular in younger clades (the “out of equilibrium” hypothesis, proposed to explain the negative values of $$\gamma$$ often observed in real phylogenies; Liow et al., 2010; Gascuel et al., 2015; Manceau et al., 2015; Missa et al., 2016; Bonnet-Lebrun et al., 2017). Taking into account such correlations between the age of clades and their branching tempo would also affect the expected loss of PD. The EDGE program (“Evolutionary Distinct and Globally Endangered”; Isaac et al., 2007) encourages conservation priorities aiming at preserving most evolutionary history within the tree of life, by proposing a ranking of species based on combined criteria of evolutionary distinctiveness and extinction risk. Although our approach is not species-based but clade-based, it also investigates the preservation of evolutionary history based on principles linked to species evolutionary distinctiveness (related to the depths of subtrees, which depend on $$\alpha$$) and to the distribution of extinction risks in the tree (which depends on $$\eta$$). Accordingly, to conserve most evolutionary history and evolutionary potential for further diversification and/or survival, priority could be given to clades that would undergo higher loss of PD in the face of species extinctions, that is, clades in the danger zone ($$\eta>1$$ and either $$\beta<-1$$ or $$\alpha<0)$$, and although not shown but only discussed herein, with $$\gamma<0$$ (decreasing branching tempo; Pybus and Harvey, 2000). Beyond Losses of PD As we have seen earlier, the trait value of a species initially interpreted as its relative abundance can also be seen as a probability of sampling it. So an additional tool provided by the parameter $$\eta$$ is a sampling distribution where sampling is not independent of the phylogeny. In particular, our results could be interpreted in the light of rarefaction experiments (Nipperess and Matsen, 2013), which study the way phylogenetic patterns in a metacommunity change as sampling decreases. Previous studies already pointed out strong impacts of nonrandom taxon sampling on the macroevolutionary patterns that we observe (e.g., Cusimano and Renner, 2010). Our results provide insights on the effects of nonrandom sampling on PD and phylogenetic tree topology. They reveal how, when the rarest species are not sampled, the discrepancy between observed and real PD depends on tree balance and abundance-richness relation (this discrepancy being larger in particular in the danger zone; Fig. 6E); and how the discrepancy between observed and real tree shape depends on $$\eta$$ (real trees being more unbalanced if $$\eta<1$$, and diverging from Yule trees towards more balance or more unbalance if $$\eta>1$$; Fig. 8). These effects of incomplete sampling on macroevolutionary patterns should be particularly important to understand biodiversity patterns in bacterial and archeal phyla, which remain poorly known in particular because they likely harbor rare species having high chances to remain unnoticed. Conclusion This new stochastic model of phylogenetic trees spans a large range of binary trees endowed with node rankings and species abundances/range sizes/extinction risks, based on three parameters only and interpolating other well-known models. We showed that ranked tree shapes, nonrandom extinctions and the interactions thereof, may have a strong impact on the loss of PD in the face of species extinctions, potentially reversing some expected patterns of variation in PD. The simplicity of the model allows one to infer the parameters on empirical phylogenies. Applying our inference procedure on bird family phylogenies we found that, in this data set, the parameters fall within a narrow range of the parameter space; and that the inferred values make the PD of these trees very sensitive to species extinctions. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.mv980. Acknowledgments The authors are very grateful to Mike Steel, Ana S.L. Rodrigues, and Hélène Morlon for their feedbacks on an earlier version of this manuscript. They are indebted to Arne Mooers and to Walter Jetz for providing the data on bird range sizes. They wish to thank the Associate Editor and the referees for their careful reviews. They thank the Center for Interdisciplinary Research in Biology (Collège de France, CNRS) for funding. References Aldous D. 1996. Probability distributions on cladograms. In: Aldous, D., Pemantle, R., editors. Random discrete structures . New York: Springer. p. 1– 18. Google Scholar CrossRef Search ADS   Aldous D. 2001. Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Stat. Sci.  16: 23– 34. Google Scholar CrossRef Search ADS   Alexander H.K., Lambert A., Stadler T. 2015. Quantifying age-dependent extinction from species phylogenies. Syst. Biol.  65: 35– 50. Google Scholar CrossRef Search ADS PubMed  Alfaro M., Santini F., Brock C., Alamillo H., Dornburg A., Rabosky D., Carnevale G., Harmon L. 2009. Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates. Proc. Natl. Acad. Sci. USA  106: 13410– 13414. Google Scholar CrossRef Search ADS   Arbetman, M.P., Gleiser G., Morales C.L., Williams P., Aizen M.A. 2017. Global decline of bumblebees is phylogenetically structured and inversely related to species range size and pathogen incidence. Proc. Roy. Soc. Lond. B  284: 20170204. Google Scholar CrossRef Search ADS   Baillie, J.E.M., Hilton-Taylor C., Stuart S.N. 2004. A global species assessment. Gland, Switzerland: IUCN. Google Scholar CrossRef Search ADS   Barnosky, A.D., Matzke N., Tomiya S., Wogan G.O.U., Swartz B., Quental T.B., Marshall C., McGuire J.L., Lindsey E.L., Maguire K.C., Mersey B., Ferrer E.A. 2011. Has the Earth’s sixth mass extinction already arrived? Nature  471: 51– 57. Google Scholar CrossRef Search ADS PubMed  Barraclough, T.G. 2010. Evolving entities: towards a unified framework for understanding diversity at the species and higher levels. Philos. Trans. R. Soc. Lond. B  365: 1801– 1813. Google Scholar CrossRef Search ADS   Bennett, P.M., Owens I.P.F. 1997. Variation in extinction risk among birds: chance or evolutionary predisposition? Proc. R. Soc. Lond. B.  264: 401– 408. Google Scholar CrossRef Search ADS   Bertoin, J. 2002. Self-similar fragmentations. Ann. I.H. Poincaré  38: 319– 340. Google Scholar CrossRef Search ADS   Bertoin, J. 2006. Random fragmentation and coagulation processes. Cambridge: Cambridge University Press. Google Scholar CrossRef Search ADS   Bielby, J., Cunningham A.A., Purvis A. 2006. Taxonomic selectivity in amphibians: ignorance, geography or biology? Anim.  Conserv.  9: 135– 143. Blum, M., François O. 2006. Which random processes describe the tree of life? A large-scale study of phylogenetic tree imbalance. Syst. Biol.  55: 685– 691. Bokma, F. 2003. Testing for equal rates of cladogenesis in diverse taxa. Evolution  57: 2469– 2474. Google Scholar CrossRef Search ADS PubMed  Bonnet-Lebrun, A.-S., Manica A., Eriksson A., Rodrigues A.S. 2017. Empirical phylogenies and species abundance distributions are consistent with preequilibrium dynamics of neutral community models with gene flow. Evolution  71: 1149– 1163. Google Scholar CrossRef Search ADS PubMed  Bortolussi N., Durand E., Blum M., François O. 2006. apTreeshape: statistical analysis of phylogenetic tree shape. Bioinformatics  22: 363– 364. Google Scholar CrossRef Search ADS PubMed  Cardillo, M., Mace G.M., Gittleman J.L., Purvis A. 2006. Latent extinction risk and the future battlegrounds of mammal conservation. Proc. Roy. Soc. Lond. B  103: 4157– 4161. Cardillo M., Orme C.D.L., Owens I.P.F. 2005. Testing for latitudinal bias in diversification rates : an example using New World birds. Ecology  86: 2278– 2287. Google Scholar CrossRef Search ADS   Colwell R.K., Lees D.C. 2000. The mid-domain effect: geometric species richness. Trends Ecol. Evol.  15: 70– 76. Google Scholar CrossRef Search ADS PubMed  Cusimano N., Renner S.S. 2010. Slowdowns in diversification rates from real phylogenies may not be real. Syst. Biol.  59: 458– 464. Google Scholar CrossRef Search ADS PubMed  Cusimano N., Stadler T., Renner S.S. 2012. A new method for handling missing species in diversification analysis applicable to randomly or nonrandomly sampled phylogenies. Syst. Biol.  61: 785– 792. Google Scholar CrossRef Search ADS PubMed  Davies T.J., Smith G.F., Bellstedt D.U., Boatwright J.S., Bytebier B., Cowling R.M., Forest F., Harmon L.J., Muasya A.M., Schrire B.D., Steenkamp Y., van der Bank M., Savolainen V. 2011. Extinction risk and diversification are linked in a plant biodiversity hotspot. PLoS Biol.  9: 1– 9. Google Scholar CrossRef Search ADS   Davies T.J., Yessoufou K. 2013. Revisiting the impacts of non-random extinction on the tree-of-life. Biol. Lett.  9: 20130343. Google Scholar CrossRef Search ADS PubMed  Doran N.A., Arnold A.J., Parker W.C., Huffer F.W. 2006. Is extinction age dependent? Palaios  21: 571– 579. Google Scholar CrossRef Search ADS   Ewens W.J. 2012. Mathematical population genetics 1: theoretical introduction. New York: Springer Science & Business Media. Google Scholar CrossRef Search ADS   Faith D.P. 1992. Conservation evaluation and phylogenetic diversity. Biol. Conserv.  61: 1– 10. Google Scholar CrossRef Search ADS   Faller B., Pardi F., Steel M. 2008. Distribution of phylogenetic diversity under random extinction. J. Theor. Biol.  251: 286– 296. Google Scholar CrossRef Search ADS PubMed  Feng S. 2010. The Poisson–Dirichlet distribution and related topics: models and asymptotic behaviors. Heidelberg: Springer Science & Business Media. Google Scholar CrossRef Search ADS   Forest F. 2009. Calibrating the tree of life: fossils, molecules and evolutionary timescales. Ann. Bot.  104: 789– 794. Google Scholar CrossRef Search ADS PubMed  Frishkoff L.O., Karp D.S., M’Gonigle L.K., Mendenhall C.D., Zook J., Kremen C., Hadly E.A., Daily G.C. 2014. Loss of avian phylogenetic diversity in neotropical agricultural systems. Science  345: 1343– 1346. Google Scholar CrossRef Search ADS PubMed  Fritz S.A., Purvis A. 2010. Selectivity in mammalian extinction risk and threat types: a new measure of phylogenetic signal strength in binary traits. Conserv. Biol.  24: 1042– 1051. Google Scholar CrossRef Search ADS PubMed  Gascuel F., Ferrière R., Aguilée R., Lambert A. 2015. How ecology and landscape dynamics shape phylogenetic trees. Syst. Biol.  64: 590– 607. Google Scholar CrossRef Search ADS PubMed  Gaston K.J., Blackburn T.M. 1995. Birds, body size and the threat of extinction. Philos. Trans. R. Soc. Lond. B  347: 205– 212. Google Scholar CrossRef Search ADS   Glavin T. 2007. The sixth extinction: journeys among the lost and left behind. New York: Thomas Dunne Books. Guyer C., Slowinski J.B. 1991. Comparisons of observed phylogenetic topologies with null expectations among three monophyletic lineages. Evolution  45: 340– 350. Google Scholar CrossRef Search ADS PubMed  Guyer C., Slowinski J.B. 1993. Adaptive radiation and the topology of large phylogenies. Evolution  47: 253– 263. Google Scholar CrossRef Search ADS PubMed  Hagen O., Hartmann K., Steel M., Stadler T. 2015. Age-dependent speciation can explain the shape of empirical phylogenies. Syst. Biol.  64: 432– 440. Google Scholar CrossRef Search ADS PubMed  Heard S.B. 1992. Patterns in tree balance among cladistic, phenetic, and randomly generated phylogenetic trees. Evolution  46: 1818– 1826. Google Scholar CrossRef Search ADS PubMed  Heard S.B., Mooers A.O. 2000. Phylogenetically patterned speciation rates and extinction risks change the loss of evolutionary history during extinctions. Proc. R. Soc. Lond. B.  267: 613– 620. Google Scholar CrossRef Search ADS   Heath T.A., Hedtke S.M., Hillis D.M. 2008. Taxon sampling and the accuracy of phylogenetic analyses. J. Syst. Evol.  46: 239– 257. Hubbell S. 2001. The unified neutral theory of biodiversity and biogeography. Princeton, NJ: Princeton University Press. Google Scholar CrossRef Search ADS   Hughes A.L. 1999. Differential human impact on the survival of genetically distinct avian lineages. Bird Conserv. Int.  9: 147– 154. Google Scholar CrossRef Search ADS   Isaac N.J.B., Turvey S.T., Collen B., Waterman C., Baillie J.E.M 2007. Mammals on the EDGE: conservation priorities based on threat and phylogeny. PLoS One  2: e296. Google Scholar CrossRef Search ADS PubMed  IUCN. 2012. IUCN red list categories and criteria: version 3.1. 2nd ed. Gland, Switzerland: IUCN. Jetz W., Thomas G., Joy J., Hartmann K., Mooers A. 2012. The global diversity of birds in space and time. Nature  491: 444– 448. Google Scholar CrossRef Search ADS PubMed  Johnson C.N., Delean S., Balmford A. 2002. Phylogeny and the selectivity of extinction in Australian marsupials. Anim. Conserv.  5: 135– 142. Google Scholar CrossRef Search ADS   Johnson S.G., Narasimhan B. 2013. R package “cubature”: adaptive multivariate integration over hypercubes. Available from: https://cran.r-project.org/web/packages/cubature/index.html. Jønsson, K.A., Fabre P.-h., Fritz S.A., Etienne R.S., Ricklefs R.E., Jørgensen T.B. 2012. Ecological and evolutionary determinants for the adaptive radiation of the Madagascan vangas. Proc. Natl. Acad. Sci. USA  109: 6620– 6625. Google Scholar CrossRef Search ADS   Kembel S.W., Ackerly D.D., Blomberg S.P., Cornwell W.K., Cowan P.D., Helmus M.R., Morlon H., Webb C.O. 2014. R package “picante”: R tools for integrating phylogenies and ecology. Available from: https://cran.r-project.org/web/packages/picante/index.html. Kingman J.F.C. 1982. The coalescent. Stoch. Proc. Appl.  13: 235– 248. Google Scholar CrossRef Search ADS   Kirkpatrick M., Slatkin M. 1993. Searching for evolutionary patterns in the shape of a phylogenetic tree. Evolution  47: 1171– 1181. Google Scholar CrossRef Search ADS PubMed  Kozak K.H., Weisrock D.W., Larson A. 2006. Rapid lineage accumulation in a non-adaptive radiation: phylogenetic analysis of diversification rates in eastern North American woodland salamanders (Plethodontidae: Plethodon). Proc. R. Soc. Lond. B.  273: 539– 546. Google Scholar CrossRef Search ADS   Kumar S. 2005. Molecular clocks: four decades of evolution. Nat. Rev. Genet.  6: 654– 662. Google Scholar CrossRef Search ADS PubMed  Lambert A., Stadler T. 2013. Birth-death models and coalescent point processes: the shape and probability of reconstructed phylogenies. Theor. Popul. Biol.  90: 113– 128. Google Scholar CrossRef Search ADS PubMed  Lambert A., Steel M. 2013. Predicting the loss of phylogenetic diversity under non-stationary diversification models. J. Theor. Biol.  337: 111– 124. Google Scholar CrossRef Search ADS PubMed  Lambert A.et al.   2017. Probabilistic models for the (sub) tree (s) of life. Braz. J. Probab. Stat.  31: 415– 475. Google Scholar CrossRef Search ADS   Leakey R.E., Lewin R. 1995. The sixth extinction: patterns of life and the future of humankind. New York: Doubleday. Lee T.M., Jetz W. 2011. Unravelling the structure of species extinction risk for predictive conservation science. Proc. Roy. Soc. Lond. B.  278 1710: 1329– 1338 Google Scholar CrossRef Search ADS   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: 646– 659. Google Scholar CrossRef Search ADS PubMed  Lovette I.J., Bermingham E. 1999. Explosive speciation in the New World Dendroica warblers. Proc. R. Soc. Lond. B.  266: 1629– 1636. Google Scholar CrossRef Search ADS   Lozano F.D., Schwartz M.W. 2005. Patterns of rarity and taxonomic group size in plants. Biol. Conserv.  126: 146– 154. Google Scholar CrossRef Search ADS   MacArthur R., Wilson E. 1967. The theory of island biogeography. Princeton, NJ: Princeton University Press. Google Scholar CrossRef Search ADS   MacArthur R.H. 1957. On the relative abundance of bird species. Proc. Natl. Acad. Sci. USA  43: 293– 295. Google Scholar CrossRef Search ADS   Maddison W.P., Midford P.E., Otto S.P. 2007. Estimating a binary character’s effect on speciation and extinction. Syst. Biol.  56: 701– 710. Google Scholar CrossRef Search ADS PubMed  Magallon S., Sanderson M.J. 2001. Absolute diversification rates in Angiosperm clades. Evolution  55: 1762– 1780. Google Scholar CrossRef Search ADS PubMed  Magnuson-Ford K., Ingram T., Redding D.W., Mooers A.Ø. 2009. Rockfish (sebastes) that are evolutionarily isolated are also large, morphologically distinctive and vulnerable to overfishing. Biol. Conserv.  142: 1787– 1796. Google Scholar CrossRef Search ADS   Manceau M., Lambert A., Morlon H. 2015. Phylogenies support out-of-equilibrium models of biodiversity. Ecol. Lett.  18: 347– 356. Google Scholar CrossRef Search ADS PubMed  McKinney M.L. 1997. Extinction vulnerability and selectivity: combining ecological and paleontological views. Annu. Rev. Ecol. Syst.  28: 495– 516. Google Scholar CrossRef Search ADS   McPeek M.A. 2008. The ecological dynamics of clade diversification and community assembly. Am. Nat.  172: E270– E284. Google Scholar CrossRef Search ADS PubMed  McPeek M.A., Brown J.M. 2007. Clade age and not diversification rate explains species richness among animal taxa. Am. Nat.  169: E97– 106. Google Scholar CrossRef Search ADS PubMed  Missa O., Dytham C., Morlon H. 2016. Understanding how biodiversity unfolds through time under neutral theory. Philos. Trans. R. Soc. Lond. B  371: 1– 12. Google Scholar CrossRef Search ADS   Moen D., Morlon H. 2014. Why does diversification slow down? Trends Ecol. Evol.  29: 190– 197. Google Scholar CrossRef Search ADS PubMed  Mooers A., Gascuel O., Stadler T., Li H., Steel M. 2012. Branch lengths on birth-death trees and the expected loss of phylogenetic diversity. Syst. Biol.  61: 195– 203. Google Scholar CrossRef Search ADS PubMed  Mooers A., Heard S. 1997. Inferring evolutionary process from phylogenetic tree shape. Q. Rev. Biol.  72: 31– 54. Google Scholar CrossRef Search ADS   Mooers A.O. 1995. Tree balance and tree completeness. Evolution  49: 379– 384. Google Scholar CrossRef Search ADS PubMed  Mooers A. Ø., Goring S.J., Turvey S.T., Kuhn T.S. 2009. Holocene extinctions and the loss of feature diversity. Turvey S. editor. Holocene extinctions.  Oxford: Oxford University Press. p. 279– 338 Google Scholar CrossRef Search ADS   Morlon H., Parsons T.L., Plotkin J.B. 2011a. Reconciling molecular phylogenies with the fossil record. Proc. Natl. Acad. Sci. USA  108: 16327– 16332. Google Scholar CrossRef Search ADS   Morlon H., Schwilk D.W., Bryant J.A., Marquet P.A., Rebelo A.G., Tauss C., Bohannan B.J., Green J.L. 2011b. Spatial patterns of phylogenetic diversity. Ecol. Lett.  14: 141– 149. Google Scholar CrossRef Search ADS   Nee S. 2006. Birth-death models in macroevolution. Annu. Rev. Ecol. Evol. Syst.  37: 1– 17. Google Scholar CrossRef Search ADS   Nee S., May R.M. 1997. Extinction and the loss of evolutionary history. Science  278: 692– 694. Google Scholar CrossRef Search ADS PubMed  Nee S., Mooers A.O., Harvey P.H. 1992. Tempo and mode of evolution revealed from molecular phylogenies. Proc. Natl. Acad. Sci. USA  89: 8322– 8326. Google Scholar CrossRef Search ADS   Nipperess D.A., Matsen F.A. 2013. The mean and variance of phylogenetic diversity under rarefaction. Methods Ecol. Evol.  4: 566– 572. Google Scholar CrossRef Search ADS PubMed  Paradis E., Claude J., Strimmer K. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics  20: 289– 290. Google Scholar CrossRef Search ADS PubMed  Parhar R.K., Mooers A.Ø. 2011. Phylogenetically clustered extinction risks do not substantially prune the tree of life. PLoS One  6: e23528. Google Scholar CrossRef Search ADS PubMed  Pearson P.N. 1992. Survivorship analysis of fossil taxa when real-time extinction rates vary: the paleogene planktonic foraminifera. Paleobiology  18: 115– 131. Google Scholar CrossRef Search ADS   Phillimore A.B., Price T.D. 2008. Density-dependent cladogenesis in birds. PLoS Biol . 6: e71. Google Scholar CrossRef Search ADS PubMed  Pigot A.L., Phillimore A.B., Owens I.P.F., Orme C.D.L. 2010. The shape and temporal dynamics of phylogenetic trees arising from geographic speciation. Syst. Biol.  59: 660– 673. Google Scholar CrossRef Search ADS PubMed  Prado P.I., Miranda M.D., Chalom A. 2015. R package “sads”: maximum likelihood models for species abundance distributions. Available form: http://search.r-project.org/library/sads/html/fitsad.html. Pulquério, M.J.F., Nichols R.A. 2007. Dates from the molecular clock: how wrong can we be? Trends Ecol. Evol.  22: 180– 184. Google Scholar CrossRef Search ADS PubMed  Purvis A. 1996. Using interspecies phylogenies to test macroevolutionary hypotheses. In: Harvey, P., Leigh Brown, A., Maynard Smith, J., Nee, S., editors. New uses for new phylogenies . Oxford: Oxford University Press. p. 153– 168. Purvis A. 2008. Phylogenetic approaches to the study of extinction. Annu. Rev. Ecol. Evol. Syst.  39: 301– 319. Google Scholar CrossRef Search ADS   Purvis A., Agapow P.-M., Gittleman J.L., Mace G.M. 2000a. Nonrandom extinction and the loss of evolutionary history. Science  288: 328– 330. Google Scholar CrossRef Search ADS   Purvis A., Gittleman J.L., Cowlishaw G., Mace G.M. 2000b. Predicting extinction risk in declining species. Proc. Roy. Soc. Lond. B  267: 1947– 1952. Google Scholar CrossRef Search ADS   Purvis A., Jones K.E., Mace G.M. 2000c. Extinction. BioEssays  22: 1123– 1133. Google Scholar CrossRef Search ADS   Pybus O.G., Harvey P.H. 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc. R. Soc. Lond. B.  267: 2267– 2272. Google Scholar CrossRef Search ADS   R Development Core Team. 2012. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Rabosky D.L. 2009. Ecological limits and diversification rate: alternative paradigms to explain the variation in species richness among clades and regions. Ecol. Lett.  12: 735– 743. Google Scholar CrossRef Search ADS PubMed  Rabosky D.L., Lovette I.J. 2008a. Density-dependent diversification in North American wood warblers. Proc. R. Soc. Lond. B.  275: 2363– 2371. Google Scholar CrossRef Search ADS   Rabosky D.L., Lovette I.J. 2008b. Explosive evolutionary radiations: decreasing speciation or increasing extinction through time? Evolution  62: 1866– 1875. Google Scholar CrossRef Search ADS   Rabosky D.L., Slater G.J., Alfaro M.E. 2012. Clade age and species richness are decoupled across the eukaryotic tree of life. PLoS Biol . 10: e1001381. Google Scholar CrossRef Search ADS PubMed  Raup D.M., Gould S.J., Schopf T.J.M., Simberloff D.S. 1973. Stochastic models of phylogeny and the evolution of diversity. J. Geol.  81: 525– 542. Google Scholar CrossRef Search ADS   Redding D.W., Mooers A.Ø. 2006. Incorporating evolutionary measures into conservation prioritization. Conserv. Biol.  20: 1670– 1678. Google Scholar CrossRef Search ADS PubMed  Ricklefs R. 2009. Speciation, extinction and diversity. In: Butlin, R., Bridle, J., Schluter, D., editors. Speciation and patterns of diversity . Cambridge: Cambridge University Press. p. 257– 277. Google Scholar CrossRef Search ADS   Ricklefs R.E. 2006. Global variation in the diversification rate of passerine birds. Ecology  87: 2468– 78. Google Scholar CrossRef Search ADS PubMed  Ricklefs R.E. 2007. History and diversity: explorations at the intersection of ecology and evolution. Am. Nat.  170: S56– S70. Google Scholar CrossRef Search ADS PubMed  Rüber L., Zardoya R. 2005. Rapid cladogenesis in marine fishes revisited. Evolution  59: 1119– 1127. Google Scholar CrossRef Search ADS PubMed  Russell G.J., Brooks T.M., McKinney M.M., Anderson C.G. 1998. Present and future taxonomic selectivity in bird and mammal extinctions. Conserv. Biol.  12: 1365– 1376. Google Scholar CrossRef Search ADS   Sainudiin R., Véber A. 2016. A beta-splitting model for evolutionary trees. R. Soc. Open Sci.  3: 160016. Google Scholar CrossRef Search ADS PubMed  Sánchez-Reyes L.L., Morlon H., Magallón S. 2016. Uncovering higher-taxon diversification dynamics from clade age and species-richness data. Syst. Biol.  66: 367– 378. Schwartz M.W., Simberloff D. 2001. Taxon size predicts rates of rarity in vascular plants. Ecol. Lett.  4 5: 464– 469. Google Scholar CrossRef Search ADS   Schwartz R.S., Mueller R.L. 2010. Branch length estimation and divergence dating: estimates of error in Bayesian and maximum likelihood frameworks. BMC Evol. Biol.  10: 1– 21. Google Scholar CrossRef Search ADS PubMed  Seehausen O. 2006. African cichlid fish: a model system in adaptive radiation research. Proc. R. Soc. Lond. B.  273: 1987– 1998. Google Scholar CrossRef Search ADS   Slowinski J., Guyer C. 1993. Testing whether certain traits have caused amplified diversification—an improved method based on a model of random speciation and extinction. Am. Nat.  142: 1019– 1024. Google Scholar CrossRef Search ADS PubMed  Stadler T. 2013. Recovering speciation and extinction dynamics based on phylogenies. J. Evol. Biol.  26: 1203– 1219. Google Scholar CrossRef Search ADS PubMed  Stadler T., Rabosky D.L., Ricklefs R.E., Bokma F. 2014. On age and species richness of higher taxa. Am. Nat.  184: 447– 455. Google Scholar CrossRef Search ADS PubMed  Szabo J.K., Khwaja N., Garnett S.T., Butchart S.H. 2012. Global patterns and drivers of avian extinctions at the species and subspecies level. PLoS One  7: e47080. Google Scholar CrossRef Search ADS PubMed  Vamosi J.C., Wilson J.R. 2008. Nonrandom extinction leads to elevated loss of angiosperm evolutionary history. Ecol. Lett.  11: 1047– 1053. Google Scholar CrossRef Search ADS PubMed  Van Valen L. 1976. Ecological species, multispecies, and oaks. Taxon  25: 233– 239. Google Scholar CrossRef Search ADS   Vazquez D.P., Gittleman J.L. 1998. Biodiversity conservation: does phylogeny matter? Curr. Biol.  8: 379– 381. Google Scholar CrossRef Search ADS   Venditti C., Meade A., Pagel M. 2010. Phylogenies reveal new interpretation of speciation and the Red Queen. Nature  463: 349– 352. Google Scholar CrossRef Search ADS PubMed  Veron S., Davies T.J., Cadotte M.W., Clergeau P., Pavoine S. 2015. Predicting loss of evolutionary history: where are we? Biol. Rev.  0: 0– 0. von Euler F. 2001. Selective extinction and rapid loss of evolutionary history in the bird fauna. Proc. R. Soc. Lond. B.  268: 127– 130. Google Scholar CrossRef Search ADS   Wake D.B., Vredenburg V.T. 2008. Are we in the midst of the sixth mass extinction? A view from the world of amphibians. Proc. Natl. Acad. Sci. USA  105: 11466– 11473. Google Scholar CrossRef Search ADS   Weir J.T. 2006. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution  60: 842– 855. Google Scholar CrossRef Search ADS PubMed  Welch J.J., Bromham L. 2005. Molecular dating when rates vary. Trends Ecol. Evol.  20: 320– 327. Google Scholar CrossRef Search ADS PubMed  Yule G. 1925. A mathematical theory of evolution, based on the conclusions of Dr. J.C. Willis, F.R.S. Philos. Trans. R. Soc. Lond. B  213: 402– 410. Google Scholar CrossRef Search ADS   Zink R., Slowinski J. 1995. Evidence from molecular systematics for decreased avian diversification in the Pleistocene Epoch. Proc. Natl. Acad. Sci. USA  92: 5832– 5835. Google Scholar CrossRef Search ADS   © 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/about_us/legal/notices)

### Journal

Systematic BiologyOxford University Press

Published: Apr 16, 2018

## 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