Odile Maliet
email: odile.maliet@orange.fr

Florian Hartig

Hélène Morlon

A model with many small shifts for estimating species-specific diversification rates

Understanding how and why diversification rates vary through time, space, and across species groups is key to understanding the emergence of today's biodiversity. Phylogenetic approaches aimed at identifying variations in diversification rates during the evolutionary history of clades have focused on exceptional shifts subtending evolutionary radiations. While such shifts have undoubtedly affected the history of life (1), identifying smaller but more frequent changes is important as well. We develop ClaDS, a new Bayesian approach for estimating branch-specific diversification rates on a phylogeny, that relies on a model with changes in diversification rates at each speciation event. We show using Monte-Carlo simulations that the approach performs well at inferring both small and large changes in diversification. Applying our approach

to bird phylogenies covering the entire avian radiation, we find that diversification rates are remarkably heterogeneous within evolutionary restricted species groups. Some groups such as Accipitridae (hawks and allies) cover almost the full range of speciation rates found across the entire bird radiation. As much as 76% of the variation in branch-specific rates across this radiation is due to intra-clade variation, suggesting that a large part of the variation in diversification rates is due to many small rather than few large shifts.

Manuscript

Several phylogenetic approaches have been developed for understanding when and on which lineages diversification rates have changed during the evolutionary history of clades (2; 1; 3; 4; 5; 6).

Most have focused on 'major' rate shifts, which is convenient methodologically; likelihoods of trees under such models have been used for some time [START_REF] Sanderson | Diversification rates in a temperate legume clade: are there so many species of Astragalus (Fabaceae)?[END_REF]. These models correspond to the idea that few rare events, such as key innovations, facilitate the invasion of new adaptive zones, with a drastic impact on diversification rates (8; 9). In these models, outside of few remarkable events, diversification rates are assumed to be homogeneous. However, while major rate shifts linked to key innovations have undoubtedly affected the history of life [START_REF] Alfaro | Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates[END_REF], they are not the only -nor necessarily the most important -source of variation in diversification rates. Shifts in diversification rates are likely quite widespread. Speciation and extinction rates may vary across lineages as a response to the particular biotic and abiotic environment experienced by each lineage [START_REF] Benton | The Red Queen and the Court Jester: species diversity and the role of biotic and abiotic factors through time[END_REF]; they may also vary as a response to traits that affect reproductive isolation such as reproduction mode [START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF] or pollination and dispersal syndromes [START_REF] Onstein | Frugivory-related traits promote speciation of tropical palms[END_REF]. Such changes in diversification rates probably occur far more frequently than key innovations, resulting in heterogeneous 1 diversification rates at much finer taxonomic scales [START_REF] Jetz | The global diversity of birds in space and time[END_REF]. Accounting for such finer scale heterogeneity is crucial if we want to obtain refined estimates of lineage-specific diversification rates and to better understand the processes subtending heterogeneity in the diversification of life. Methods of the State-Speciation-Extinction family [START_REF] Fitzjohn | Diversitree: comparative phylogenetic analyses of diversification in R[END_REF] can in principle better account for these types of heterogeneities, but they require assuming trait-dependency of rates (SI Appendix section 3.5). Nonmodel-based approaches such as the DR statistic (4) can also account for fine-scale heterogeneities, but they are rather ad hoc and generally do not perform as well as model-based approaches [START_REF] Title | Tip rates, phylogenies, and diversification: what are we estimating, and how good are the estimates? Methods in Ecology and Evolution[END_REF].

Here, we develop a new Bayesian approach (ClaDS) for estimating lineage-specific diversification rates on a phylogeny that better accounts for the diverse sources of variation in diversification rates that occur during the evolutionary history of clades. Using Monte Carlo simulations, we quantify the ability of ClaDS to faithfully recover both small and large changes in diversification rates. Finally, we apply the method to time-calibrated phylogenies for 42 bird clades to evaluate the extent to which differences in the pace of diversification across the entire avian radiation result from few large versus many small events.

A new model of diversification rate variation

We consider a birth-death diversification process, the cladogenetic diversification rate shift (ClaDS) model, where diversification rates are inherited at speciation, but with a shift (Fig. 1).

At the beginning of the process, the clade is composed of one lineage with speciation rate λ 0 and extinction rate µ 0 . At each speciation event, the two daughter lineages inherit new diversification rates (λ i1 , λ i2 ) and (µ i1 , µ i2 ) sampled from a joint probability distribution ν parameterized by the parental rates λ i and µ i . If the change in speciation and extinction rates are assumed to be independent, the λ i are sampled from a distribution ν λ , the µ i are sampled from a distribution ν µ , and ν = ν λ × ν µ . Moreover, we allow for the possibility that some extant species are missing by assuming that each extant species is observed with probability f ≤ 1. We derive the probability density of a reconstructed phylogeny under this general model and implement its computation in R (Materials and Methods & SI Appendix sections 2 to 5).

We then consider several scenarios in ClaDS where: i) ν λ is a lognormal distribution with parameters log(α * λ) and σ; the latter ensures that the relative change in rate at speciation λ i /λ is independent from the parental rate with a mean m given by α exp(σ 2 /2); σ controls how constrained daughter rates are (highly constrained for small σ values) and α controls the trend at speciation (i.e. whether daughter rates tend to be higher or lower than parental rates) ii) extinction rates are either negligible (µ i = 0 for all lineages, ClaDS0), homogeneous across all lineages in the clade (µ i = µ 0 for all lineages, ClaDS1) or vary across lineages, but with a constant turnover ε (i.e. µ i /λ i = ε for all lineages, ClaDS2). We use Monte Carlo simulations under ClaDS1 and ClaDS2 (Materials and Methods & SI Appendix section 6) to verify that our likelihood expression is correct (SI Appendix section 6, Fig. S6 to S8). Finally, we implement a Monte Carlo Markov Chain (MCMC) sampler that, given a reconstructed phylogeny, simultaneously estimates both the parameters of ClaDS (λ 0 , α, σ, and either µ 0 or ε) and the speciation rates λ i at the origin of each branch i of the phylogeny (Materials and Methods & SI Appendix section 7, see also section 8 and Fig. S9 & S10 for a test of the sampler). Branch-specific extinction rates µ i at the origin of each branch i of the phylogeny are given by µ 0 for ClaDS1 and by ε * λ i for ClaDS2. In what follows for simplicity we refer to λ i and µ i as "branch-specific rates" instead of the more accurate "rates at the origin of each branch". Under these scenarios of the ClaDS process, heterogeneity in speciation rates across lineages is determined on the one hand by a stochastic component (controlled by σ), and on the other hand by a trend component (controlled by m). When the expected daughter rate is equal to the parental rate (m = 1), the resulting trees are relatively imbalanced and tippy (SI Appendix section 1, Fig. S1 &S2): lineages that by chance have high speciation rates early in clade's history spread, leading to rates that are heterogeneous across lineages and average rates that increase through time. This sorting effect is exacerbated when the expected daughter rate is higher than the parental rate (m > 1, Fig. S1 &S2), corresponding to a 'niche-piling' scenario where diversity begets diversity [START_REF] Emerson | Species diversity can drive speciation[END_REF]. To the contrary, when the expected daughter rate is lower than the parental rate (m < 1), corresponding to a 'niche-filling' scenario where diversification gets harder as new species arise (16; 17; 18), the heterogeneity in speciation rates across lineages is reduced, and with a low enough m, the average rate is constant or even decreasing through time (Fig. S1 &S2). Importantly, ClaDS is able to produce the combination of stemmy and imbalanced tree shapes observed in nature, and under a wider set of parameter values for the scenario with constant turnover (ClaDS2) than the scenario with constant extinction rate (ClaDS1, Fig. S1 &S2).

Performance of ClaDS

We begin by testing the performance of ClaDS under frequent rate changes and in the absence of extinction (ClaDS0) (Materials and Methods). We find that the approach provides unbiased estimates of all model's parameters for large enough trees (size 200, Fig. 2); the relative change in rate at speciation m is also well estimated (Fig. 2 d). As expected, bias and variability around parameter estimates increase for smaller trees (Fig. S11 to S14).

ClaDS provides reliable estimates of branch specific speciation rates on average: while low rates tend to be slightly overestimated and large rates slightly underestimated, ClaDS can detect regions of the tree with relatively high or low rates (Fig. 3

& Fig. S15 & S16).
When considering also extinctions, focusing on the scenario with constant turnover (ClaDS2) as it generally produced tree shapes closer to those observed in nature, we found that estimates remain accurate at low levels of extinction (ε = 0.1) for both model parameters (Fig. S20) and branch-specific speciation rates (Fig. S21). At high levels of extinction (ε = 0.9), σ and, when the mean change in rate at speciation m approaches 1, branch-specific speciation rates, remain well estimated. It is not the case, however, of the turnover rate ε, α, and branch-specific speciation rates when m < 1, although accounting for extinction does improve inferences over ignoring it (Fig.

S20 & S21

). When extinction is not accounted for, estimated branch-specific speciation rates are generally lower than realized ones, but higher than realized net diversification rates (Fig. S21C &D).

If there are a small number of major rate shifts during the evolution of clades, rather than many small changes (tested here with a single rate shift, Materials and Methods), ClaDS is still able to provide reliable estimates of branch-specific rates (Fig. S17 &S19). The model is also able to detect when two branches in the tree belong to distinct speciation regimes as soon as the difference in rates between the two regimes is large enough (a two-fold increase or decrease in our simulations) and both regimes are represented by a large enough number of branches in the phylogeny (Fig. S19 left). The false detection rate associated to this test is low (Fig. S19 right).

Finally, when comparing the performance of ClaDS to that of two other popular methods for estimating branch-specific rates (the DR statistic and BAMM (4; 5)) under various simulation schemes (SI Appendix section 9), we find, overall, that ClaDS outperforms the other methods for trees simulated with both many small shifts at speciation (Fig S22 & S26 to S28) and gradual changes along branches (Fig S24 ), and that it performs as well as other methods for trees simulated with few large shifts (Fig. S23) and variations in extinction rates (Fig. S25). Importantly, ClaDS provides reliable estimates of the variance in rates under both the many small and the few large shifts scenarios (Fig. S30 & S32), while BAMM underestimates the variance in rates under the many small shifts scenario (Fig. S31). ClaDS and BAMM provide low (and similar) estimates of the variance in rates for trees simulated under constant rates (Fig. S29); in the presence of rate heterogeneity, they tend to underestimate rather than overestimate rate variance (Fig. S30 to S33), and BAMM more so than ClaDS.

Diversification across the avian radiation

When applying ClaDS to major bird clades (Materials and Methods), we found that lineagespecific speciation rates can vary by as much as 2 orders of magnitude within clades (Fig. 4e). In Accipitridae (hawks and allies) for example, speciation rates range from 0.013 to 1.2 Mya -1 , which almost covers the range found across the entire avian radiation (0.013 -5 Mya -1 ). Comparable within-clade heterogeneities occur in other clades, such as Muscicapidae & Turdidae, Tyrannidae and Parulidae (Fig. 4e, in orange). Such within-clade heterogeneities are way above heterogeneities arising from estimation error (Fig. S29). A variance partitioning of speciation rates across the bird radiation (Material and Methods) reveals that intra-clade variance accounts for 76% of the total variance. In comparison, BAMM would have estimated much less within-clade heterogeneities, with an intra-clade variance accounting for only 46% of the total variance (Fig. S34). Given our simulation results, this suggests that BAMM underestimates the intra-clade variance, and thus that many small shifts occurred during bird diversification that BAMM cannot detect.

While some clades have very heterogeneous rates, others are quite homogeneous, such as Ramphastides, Alcedinidae, Charadrii and Phasianidae (Fig. 4e, in blue). We did not find any significant relationship between the variance in rate values within a clade and the size (p = 0.49) or age (p = 0.93) of the clade, indicating that rate heterogeneity is not a mere result of time or species richness; rather, rates are pretty constrained in some old and rich clades (e.g. Phasianidae) as well as in some younger or less species-rich clades (e.g. Alcedinidae), while they can take very different values for distinct species of both old or young clades (e.g. Parulidae, Tyrannidae). The wide range of σ estimates found across bird clades (Fig. 4a), in comparison with rather tight α and m estimates (Fig. 4b &c), suggests that differences in rate heterogeneity across clades are due to the stochastic component of the model, rather than its trend component. Indeed, α ranges between 0.38 and 1.02 (with a mean of 0.71, Fig. 4b), which indicates a universal tendency for daughter rates to be smaller than ancestral ones, with a decline that is comparable in magnitude across clades. There is only one case when m is clearly above 1 (1.12 in Campephagidae); this corresponds to a case when most shifts correspond to rate declines, but the few shifts that correspond to rate increases are much bigger in magnitude.

Discussion

Models of diversification applied to phylogenies of extant taxa are increasingly used to understand the long-term evolution of biodiversity. These approaches have highlighted how much variable diversification rates can be across the tree of life, and the importance of these variations for explaining current patterns of diversity (the so-called 'diversification rate hypothesis' [START_REF] Rosenzweig | Species diversity gradients: we know more and less than we thought[END_REF]). Yet, despite recent advances in phylogenetic approaches for understanding diversification, detecting diversification rate variations and the processes underlying these variations remain a challenge spurring a heated debate (20; 21; 22; 23; 24). In this paper, we have developed ClaDS, a new model with frequent small variations in diversification rates together with a method to infer branch-specific diversification rates on a phylogeny. We have shown using simulations that ClaDS accurately estimates branch-specific rates. Finally, applying ClaDS to the bird phylogeny, we have shown that small but frequent changes have been instrumental in shaping global rate variation during the avian radiation.

One of the major advances of our model is to rely on an explicit and exact computation of the likelihood in the presence of extinction. Previous likelihood expressions under diversification models with variable rates were computed with the underlying assumption that shifts do not occur in extinct lineages (1; 3; 5), except in the case of trait-dependent models (see SI Appendix section 3.5 for further discussion); this is biologically implausible and can introduce an important bias depending on the intensity of extinction (22; 23). In ClaDS we relax this inconvenient assumption by integrating appropriate Ordinary Differential Equations (ODEs, SI Appendix section 3). This allows computing likelihoods accounting for rate shifts on extinct lineages, which has so far only been done through intense and impractical Monte Carlo simulations [START_REF] Moore | Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures[END_REF]. The ODE integration is computationally intensive, but not as much as to prevent running ClaDS on reasonably sized trees, as we illustrated on the bird phylogenies. Despite this significant improvement, our simulations show that estimating extinction remains difficult, in line with the well-known difficulty of estimating extinction from phylogenies of only extant taxa [START_REF] Rabosky | Extinction rates should not be estimated from molecular phylogenies[END_REF]. This is true even when simulations and inferences are performed under simple models with constant extinction or turnover rate. Despite difficulties in estimating extinction rates, properly accounting for extinctions in the likelihood computation is satisfying on a biological and theoretical standpoint, and, as we have shown, improves the estimation of both model parameters and branch specific speciation rates.

Another advantage of ClaDS is to avoid using model selection to select the number and location of rate shifts, by assuming that shifts happen at each speciation event. In the frequently-used MEDUSA method (1), stepwise AIC is used to perform this selection, with associated statistical limitations [START_REF] May | How well can we detect lineage-specific diversification-rate shifts? A simulation study of sequential AIC methods[END_REF]. In the approach of Morlon et al. [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF], likelihood ratio tests are performed to select the number of shifts, but the location of these shifts needs to be fixed a priori. Finally, in the popular Bayesian analysis of macroevolutionary mixtures (BAMM, 5), reversible jump mcmc is used, with a prior on the number and location of shifts that may influence the results (22; 26). ClaDS avoids these limitations, while still performing well in the presence of rare rate shifts with large effects.

Maybe more importantly than these technical aspects, ClaDS represents a view of evolution distinct from that of previous models: existing models focus on a small number of discrete diversification shift events spread across the tree, an idea that fits well with the concept of key innovations driving major diversification shifts (3; 1; 5); to the contrary, ClaDS allows for frequent variations linked, for example, to changes in environmental conditions or associations with continuously evolving heritable traits. Accordingly, ClaDS does not aim at identifying specific nodes in a phylogeny subtending major diversification rate shifts. Rather, it assumes that rate shifts happen at each speciation event and focuses on estimating branch-specific diversification rates. In nature, both many shifts with small effects and few shifts with large effects are likely to occur, and so it is reassuring to see that ClaDS can properly estimate branch specific rates under these two evolutionary processes.

Accurately estimating branch specific diversification rates is a critical step for understanding the processes that lead some species groups to diversify faster than others. For example, species' traits can modulate their propensity to diversify, and tests based on assessing the correlation between trait values at a phylogenies' tips and metrics capturing the diversification rate of the corresponding lineages ('tip-rate correlations' tests) have been developed to detect such effects [START_REF] Freckleton | Relating traits to diversification: a simple test[END_REF]. These types of tests have regained interest lately (see e.g. STRAPP [START_REF] Rabosky | A robust semi-parametric test for detecting trait-dependent diversification[END_REF], FiSSE [START_REF] Rabosky | A simple nonparametric test for the effects of a binary character on lineage diversification rates[END_REF], ES-sim (29), pNoTO (30; 31)), as an alternative or complement to state-dependent speciation-extinction (SSE) methods that jointly model diversification dynamics and trait evolution (32; 13). However, current metrics of species-level diversification rates have limitations. Some of them are derived from BAMM [START_REF] Rabosky | Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees[END_REF] and thus reflect a limited set of diversification rate regimes rather than lineage-specific rates per se.

Others are summary statistics describing phylogenetic branching patterns, such as the "node density" [START_REF] Freckleton | Relating traits to diversification: a simple test[END_REF], the "equal split" [START_REF] Redding | Incorporating evolutionary measures into conservation prioritization[END_REF], or the "diversification rate" (4) statistics; they are not rigorously derived from speciation-extinction models, and they generally perform worse than model-based approaches ( 14) (SI Fig. S22 to S27). ClaDS provides tip level estimates of diversification rates that should help identifying the specific features of a species that make it more or less prone to diversify. In the future, we could imagine a hybrid between SSE and ClaDS that would account for both traitdependent diversification and residual rate variation not accounted for by the trait, in the spirit of hidden states models (HiSSE [START_REF] Beaulieu | Detecting hidden diversification shifts in models of traitdependent speciation and extinction[END_REF], MSBD ( 6)). This could for example be done by imputing in ClaDS specific trend parameters α corresponding to trait shifts.

Changes in biotic and abiotic conditions can also modulate the tempo of diversification, leading diversification to be faster during some time periods than others. ClaDS accommodates temporal trends in rate variation, without the need to specify a specific form for this variation a priori as in time-dependent diversification models (35; 16; 3), and with more flexibility than models where a discrete rate shift at a given time point affects the whole clade [START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF]. In the future, the trend parameter α could depend on measured environmental variables; this would allow directly testing for an effect of these environmental variables on diversification, as in environment-dependent diversification models (37; 38), while accounting for residual rate variation.

Our ClaDS analysis of the avian radiation reveals a series of compelling results. First, and even though these estimates need to be taken with caution, we find significant (non-zero) turnover rates.

Second, we find a pervasive pattern of declines in speciation rates over time congruent with previous studies (16; 17; 18). Third, we find a remarkable heterogeneity in speciation rates, with per-lineage rates that vary by two orders of magnitude (0.01 -5 Mya -1 ), peaking around 0.15 Mya -1 . Fourth, we find that variability in speciation rates can be as high within than between clades, suggesting that rate variation may be much more widespread than currently thought and implemented in existing models. Finally, we highlight a remarkable difference across clades in terms of how constrained their diversification rates are, with plovers and allies on one extreme, and hawks and allies on the other extreme of a continuum between rates that vary less than 2 fold to more than 80 folds (Fig. 4e,f). These differences in how constrained diversification rates are striking and remain to be explained: these could be linked to differences in genetic architecture, developmental constraints, or biogeographies, for example. Together, our results refute the idea that speciation may be clock-like [START_REF] Hedges | Tree of life reveals clock-like speciation and diversification[END_REF] and emphasize the need to consider diversification models that embrace the pervasive heterogeneity of the evolutionary process. Further, they promise a bright future for approaches, such as ours, that relax the speciation clock similarly to the way the molecular clock has been relaxed (40; 41; 42): similar to molecular rates, diversification rates vary according to many small shifts.

Material and Methods

Likelihood, simulation and Bayesian implementation of ClaDS

Likelihood We derived the probability density of observing a reconstructed phylogeny with branches delimited by the times (t i , s i ) i∈ 1,N and speciation and extinction rates λ i and µ i at time t i (i.e. at the origin of each branch) under the cladogenetic diversification rate shift model (SI Appendix section 2 to 4). We note Θ the parameters of the new rate distribution ν. The probability density can be derived from three main probability functions: Φ Θ,λ,µ (t), the probability that a lineage alive at time t has speciation and extinction rates λ and µ and no descendant in the reconstructed phylogeny; χ Θ,λ,µ (t), the probability that a lineage alive at time t has speciation and extinction rates λ and µ and exactly one descendant species sampled in the reconstructed phylogeny; and ξ Θ,λ,µ (t, s, λ 1 , λ 2 , µ 1 , µ 2 ), the probability that a lineage alive at time t has speciation and extinction rates λ and µ and gives birth at time s to two daughter lineages that respectively have speciation rates λ 1 and λ 2 and extinction rates µ 1 and µ 2 . We obtained ordinary differential equations (ODEs) to solve for Φ, χ and ξ by considering the different events that can happen during a short time interval ∆ t and making ∆ t tend to 0 (SI Appendix section 3.1 to 3.3). Under a pure-birth model and for a completely sampled phylogeny, the ODEs can be solved analytically (SI Appendix section 4). In the presence of extinction and/or if there are missing taxa in the phylogeny, Φ, χ and ξ are computed by integrating the ODEs numerically, which is more computationally intensive (SI Appendix section 5).

Simulation We implemented a simulation algorithm of ClaDS in the R-package RPANDA (43, function sim_ClaDS) (SI Appendix section 1). In this implementation, the speciation rates of daughter lineages are drawn independently from a distribution ν λ . Their extinction rates are either drawn from a distribution ν µ , given by µ 0 (constant extinction rate scenario, ClaDS1), or given by ε * λ s i,1 and ε * λ s i,2 (constant turnover scenario, ClaDS2). ν λ and ν µ can be normal, log-normal, or uniform distributions. The simulations are continued until a stopping criterion is met, either a fixed time or a fixed number of species. In addition, sim_ClaDS takes as one of its arguments a parameter p controlling the probability that a shift happens at each speciation event (the default value p = 1 corresponds to the model investigated here), and a parameter n, controlling a maximum number of shifts (the default value n = INF corresponds to the model investigated here; if n takes a finite value, then p switches to 0 as soon as n switches have occurred).

Bayesian implementation We implemented a Bayesian inference approach for fitting ClaDS to reconstructed phylogenies in the R-package RPANDA (43, function fit_ClaDS) (SI Appendix 1 section 7). In order to fit ClaDS0 (no extinction), we use a Metropolis within Gibbs MCMC (Monte Carlo Markov Chain) sampler with a Bactrian proposal [START_REF] Yang | Searching for efficient Markov chain Monte Carlo proposal kernels[END_REF], and convergence is monitored by running three MCMC chains in parallel and computing Gelman statistics [START_REF] Gelman | Bayesian data analysis[END_REF]. In order to fit ClaDS1 and ClaDS2 (i.e. in the presence of extinction), and/or if there are missing taxa in the phylogeny, we use the faster blocked Differential Evolution (DE) MCMC sampler, with sampling from the past of the chains [START_REF] Ter Braak | Differential evolution Markov chain with snooker updater and fewer chains[END_REF]. We also ran three chains. For both with and without extinctiion, we use an inverse gamma prior with shape parameter 1 and rate parameter 0.1 for σ and a flat prior for all other parameters. Each estimate was computed as the mean over the iterations and the three chains. phylogenies of size 200 until we had a good coverage of subclade new rate/size combination (from 300 to 500 phylogenies per parameter set). In such simulations, there are only two distinct rates across the tree: the background rate and the new rate. We then ran ClaDS on each simulated tree using our run_ClaDS0 function and compared the retrieved estimates of branch-specific speciation rates for each tree to their simulated values by performing linear regressions and computing relative errors. Finally, we tested whether the model is able to detect if two branches in the tree belong to the same or distinct speciation regime(s): two branches were considered to have significantly different rates (distinct regime) if the difference in the estimated speciation rates between the two branches was of constant sign on at least 95% of the MCMC chains. We assessed the significance of speciation rate differences (and the corresponding sign) for all pairs of branches in the simulated trees. Finally, we quantified the 'proper detection' rate as the proportion of pairs for which a significant difference was inferred when the two branches indeed belonged to distinct speciation regimes (i.e. one had the background speciation rate and the other one had the new rate), and the 'false detection' rate as the proportion of pairs for which a significant difference was inferred, while the two branches actually belonged to the same speciation regime (i.e. both had either the background speciation rate or the new rate).

Diversification of the avian radiation

We applied ClaDS, accounting for extinction (ClaDS2, model with constant turnover) and incomplete sampling, to bird phylogenies. We used the MCC trees from Jetz et al. ( 4) with only the species for which there was molecular data, along with the associated sampling fractions provided by the authors. Most of these are family level phylogenies, with some spawning two or a few more families. We ran the model on the 42 bird phylogenies with more than 50 species. We report the distribution of branch-specific speciation rates across the 42 clades, as well as individual distributions for each clade. We partitioned the total variance of the logarithm of the branch specific speciation rates ( i ln(λ i ) -ln(λ)

2

, where ln(λ) is the mean of the log of the speciation rates for all branches in all clades) between the intra-clade ( i ln(λ i ) -ln(λ c i )

2

, where c i is the clade to which branch i belongs and ln(λ c ) is the mean of the log of the speciation rates for all branches in clade c) and inter-clade variance ( i ln(λ c i ) -ln(λ)

2

). We also tested for a potential correlation between the variance in rates and the size (number of tips) and age (crown age) of clades using PGLS (47) (two-sided test) on the Hacket backbone phylogeny provided in Jetz et al. [START_REF] Jetz | The global diversity of birds in space and time[END_REF].

Figure 4

 4 Figure 4 Patterns of diversification across 42 bird clades. Distributions across clades of (a) σ, (b) α, (c) m = α exp(σ 2 /2), and (d) ε values estimated with ClaDS. e: Distributions of branch specific speciation rates for each specific clade (grey and colored lines) and all clades pooled together (thick black line). Red: Accipitridae; Orange: Muscicapidae & Turdidae, Tyrannidae and Parulidae; Dark blue: Charadrii; Medium blue: Ramphastides; Light blue : Alcedinidae and Phasianidae); Brown: Scolopaci; Green: Anatinae f: Exemplar phylogenies colored according to their inferred branchspecific speciation rates, in Myr -1 , and plotted on the same time scale. Top panel: the Accipitridae phylogeny subtends very variable rates that tend to decrease through time (inferred parameters: σ = 0.67, α = 0.61, m = 0.76 and ε = 0.02). Bottom panel: the Ramphastides phylogeny subtends rather homogeneous rates (σ = 0.16, α = 0.97, m = 0.98 and ε = 0.05).

Acknowledgements

The authors are very grateful to Leandro Arístide, Julien Clavel, Jonathan Drury, Carmelo Fruciano, Sophia Lambert, Eric Lewitus, Marc Manceau, Olivier Missa, Benoît Perez, Ana Catarina Silva and Guilhem Sommeria-Klein for their helpful comments on an earlier version of this manuscript. This work was supported by an AMX grant (from Ecole Polytechique) and the Labex MemoLife to OM, PROCOPE mobility grant 57134817 to FH and HM, and the European Research Council [ERC 616419-PANDA] to HM.

Data availability

The simulated phylogenies used to test the method are available at https://github.com/OdileMaliet/ClaDS/tree/ in the file named trees.zip. All the empirical data used for the analysis were obtained from Jetz et al. (2012) study, and are available on https://www.nature.com/articles/nature11631.

Code availability

The R functions used to simulate and fit the model are available in the RPANDA R-package.

All the codes used to test our method are available on the github repository https://github.com/ OdileMaliet/ClaDS.git.

Testing the performance of ClaDS

We performed intensive simulations to test the performance of ClaDS. We tested both the performance of ClaDS under data generated by this model, and its performance for data generated with a discrete speciation rate shift. In order to assess the performance of ClaDS under a large parameter set and for a variety of tree sizes, we considered primarily the pure birth model with completely sampled phylogenies. We also considered the model with extinction and/or missing taxa, but only in a limited, computationally tractable, set of simulations.

Many small rate shifts (ClaDS model) For each combination of the following parameter values, we simulated 20 pure birth trees, stopping the simulation when a target tip number of 50, 100 and 200 was reached. λ 0 was fixed at 0.1, σ was taken in {0, 0.1, 0.18, 0.26, 0.34, 0.41}, and α in {1.2, 1.1, 1, 0.95, 0.9, 0.7}. We recorded the realized speciation rate on each branch in each of these simulations. We then ran ClaDS on each simulated tree using our run_ClaDS0 function. Lastly, we compared the retrieved estimates of λ 0 , σ and α to their simulated values; we also compared the retrieved estimates of branch-specific speciation rates for each tree to their realized values by performing linear regressions and computing relative errors (ratio of estimated versus realized rates).

In order to explore the model accounting for extinction, we simulated 5 trees of size 100 under 4 scenarios with constant turnover rate (ClaDS2), and for each condition either low (ε = 0.1) or high (ε = 0.9) turnover (8 scenarios in total). We focused on the scenario with constant turnover, because this scenario produced tree shapes similar to those of empirical trees under a wider set of parameter values than the alternative scenario with constant extinction rate (Fig. S2.1 versus S2.2). Maintaining a balance where extinction is neither negligible nor driving clades to extinction is also easier under ClaDS2. The four scenarios were as follows: i) high heterogeneity and decreasing rates : λ 0 = 0.1, σ = 0.7, α = 0.7 (mean relative change m = 0.9), ii) no heterogeneity and constant rates (equivalent to constant rate birth-death trees) :

iii) Low heterogeneity and no average change in rate at speciation : λ 0 = 0.1, σ = 0.2, α = 0.98 (m = 1) iv) Low heterogeneity and decreasing rates : λ 0 = 0.1, σ = 0.2, α = 0.88 (m = 0.9). We recorded the realized speciation rate at the beginning of each branch in each of these simulations.

We then ran ClaDS on each simulated tree using our R function, both accounting (run_ClaDS) and not accounting (run_ClaDS0) for extinction, the latter to evaluate the bias resulting from not accounting for extinction when it occurs. Lastly, we compared the retrieved estimates of σ, α, m and ε for each tree to their simulated values. We did not compare the retrieved estimates of λ 0 to the simulated values, because the estimates correspond to the speciation rate at the crown while the simulated values correspond to the speciation rate at the stem. These two rates can be very different in the presence of extinction. We also compared the retrieved estimates of branch-specific speciation rates and net diversification rates (speciation minus extinction) for each tree to their realized values by performing linear regressions and computing relative errors.

Few large rate shifts We also tested the behavior of ClaDS under a 'key innovation' scenario with only a single large rate shift during the history of the clade. In order to simulate this scenario, we used our sim_ClaDS function with λ 0 (the background rate in this case) fixed at 0.1, p (the probability that a rate shift happens at each speciation event) fixed at 0.

Competing interests

The authors declare no competing interests.

Author contribution

OM, FH and HM designed the study and performed research. OM contributed new analytical tools and analysed data. OM, FH and HM wrote the paper. cartoon phylogeny simulated under ClaDS, with branches colored according to their speciation rate (red: high rate, blue: low rate). Speciation rates are inherited at speciation with a shift determined by the probability distribution ν λ (here taken to be a lognormal distribution, insert).

Red arrows indicate speciation events (and associated diversification rate shifts) that are hidden in the reconstructed phylogeny as a result of extinction.