Emmanuel Paradis
email: emmanuel.paradis@ird.fr

Pablo A Tedesco

Bernard Hugueny

Running head: VARIATION IN SPECIATION AND EXTINCTION QUANTIFYING VARIATION IN SPECIATION AND EXTINCTION RATES WITH CLADE DATA

name 'clade data', that are the ages of clades together with their numbers of species.

We develop a general approach for the statistical modeling of variation in speciation and extinction rates, including temporal variation, unknown variation, and linear and nonlinear modeling. We show how this approach can be generalized to a wide range of situations, including testing the effects of life-history traits and environmental variables on diversification rates. We report the results of an extensive simulation study to assess the performance of some statistical tests presented here as well as of the estimators of speciation and extinction rates. These latter results suggest the possibility to estimate correctly extinction rate in the absence of fossils. An example with data on fish is presented.

KEY WORDS: birth-death models, extinction, maximum likelihood, speciation, stem ages.

The study of the tempo and mode of evolution has experienced a new wave of interest from evolutionists using new mathematical and statistical tools to analyze molecular phylogenies [START_REF] Sanderson | Reconstructing shifts in diversification rates on phylogenetic trees[END_REF][START_REF] Ricklefs | Estimating diversification rates from phylogenetic information[END_REF]. Following some initial breakthrough (e.g., [START_REF] Nee | Tempo and mode of evolution revealed from molecular phylogenies[END_REF][START_REF] Nee | The reconstructed evolutionary process[END_REF], significant progress has been achieved in biologically relevant statistical modeling of diversification, such as quantifying temporal variation in diversification [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF][START_REF] Hallinan | The generalized time variable reconstructed birth-death process[END_REF] or assessing the effects of biological traits on speciation and extinction rates [START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF][START_REF] Fitzjohn | Estimating traitdependent speciation and extinction rates from incompletely resolved phylogenies[END_REF][START_REF] Fitzjohn | Quantitative traits and diversification[END_REF]. Recent advances have also been accomplished in integrating molecular and fossil data (e.g., [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF][START_REF] Didier | The reconstructed evolutionary process with the fossil record[END_REF].

Most of these recent statistical developments have focused on analyzing complete phylogenies. Incomplete phylogenies are often treated as a seperate case in order to take missing data into account [START_REF] Pybus | New inferences from tree shape: numbers of missing taxa and population growth rates[END_REF][START_REF] Fitzjohn | Estimating traitdependent speciation and extinction rates from incompletely resolved phylogenies[END_REF][START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF]. The most common form of such data is a phylogeny resolved at a high level accompanied by the number of species associated to each tip of the tree. On the other hand, the ages of clades together with the numbers of species (named here 'clade data') have been a neglected source of data in the analysis of diversification. [START_REF] Magallón | Absolute diversification rates in angiosperm clades[END_REF] provided some methods for the analysis of such data and applied them to angiosperms. They particularly developed various estimators of the (net) rate of diversification of a clade giving its age and number of species.

The relative lack of interest towards clade data may come from the fact that, for a given clade, its complete phylogeny contains more information than the pair of values 'age + number of species'. However, for a collection of clades, such data are a valuable source of information for several reasons. First, clades defined by higher-level taxa (e.g., families, orders) are clearly identified for almost all groups of living beings and their numbers of species are in many cases already known. Second, phylogenetic relationships among higher-level taxa have been much more studied than within them, so it is more straightforward to date the age of a clade rather than the divergences among its species. Third, the fossil record is generally more informative on the origin of higher-level taxa compared to species or other low-level taxa.

Fourth, it is easier to examine the impact of the species concept on the definition of clade data rather than on a phylogeny since, in the former the species concept will mostly affect the number of species while in the latter it will be often hard to infer different phylogenies under those distinct species definitions. Clade data have also some disadvantages: the inherent lack of temporal resolution within each clade makes it impossible to study the variation in diversification within them.

In the present paper, we extend the approach presented by [START_REF] Magallón | Absolute diversification rates in angiosperm clades[END_REF] and present statistical tools for the inference of diversification patterns and processes with clade data. Our approach assumes that each clade, instead of having its own speciation and extinction rates, comes from a 'statistical population of clades' so that maximum likelihood inference is straightforward. With this rationale, we show how to make inference on variation in diversification parameters among clades using different modeling tools, including testing the effects of life-history traits and environmental variables and the case where variation is a priori unknown.

We also present the results of a simulation study in order to assess the statistical performance of several tests and estimators presented in this paper, and finally we apply our approach on a data set of fish.

Statistical Modeling Approach

Throughout this paper we assume that diversification proceeds with speciation (λ) and extinction (µ) rates which are the probabilities that a species splits into two daughter-species or goes extinct during a very short time.

We denote as X t the number of species in a clade of age t where this may be either the stem age of the clade (divergence time of the clade from its sister-clade) or its crown age (time to the most common recent ancestor of the species belonging to the clade). Specifically, using equation 8 from [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF], we can write the probability that X t takes a specific integer value x:

Pr(X t = x|θ, X 0 = 1) = η t (1 -η t ) x-1 x ≥ 1, (1)
where θ is a vector of parameters specifying how speciation and extinction rates vary through time and η t is a function of these parameters. The conditioning on X 0 = 1 emphasizes that in this paper we consider stem groups. For the case of crown groups (X 0 = 2), the probabilities must be summed on all possible combinations. In most applications, stem groups are considered because the origin of a group is inferred from its relationships with its sister-group. On the other hand, deriving the crown age of a group requires to estimate the age of the most recent common ancestor of its species which is usually more complicated because it requires to sample all species in the clade. On the other, inferring stem ages requires one species from the clade and one from its sister-clade.

Various forms exist for these probabilities depending on the parameterization of θ and whether we wish to condition them on survival of the lineage until present or not. For instance, if extinction rate is zero and speciation rate is constant, then η t = e -λt . This is the [START_REF] Yule | A mathematical theory of evolution, based on the conclusions of Dr[END_REF] model. Models with a non-null extinction rate are called birth-death models [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF].

The point of conditioning on no extinction is important when analyzing data on actual groups because total extinction of these groups did not occur.

Thus the probabilities must be modified accordingly, otherwise this would result in underestimated extinction rates [START_REF] Rabosky | Exceptional among-lineage variation in diversification rates during the radiation of Australia's most diverse vertebrate clade[END_REF]).

Let us consider for the moment the simple This approach assumes that each clade is characterized by its own speciation rate. On the other extreme, one may assume that speciation rate is the same in all groups so that the observed data are independent outcomes of the same diversification process. Thus, it is possible to use maximum likelihood inference using equation 1. The likelihood function is:

i Pr(x i |λ), (2)
where Pr(x|λ) is a simplified notation of equation 1. We may expect less bias in the estimates from this approach, but also the possibility to test hypotheses based on fitting alternative models.

The assumption of equal speciation rates among clades is, certainly in most cases, unrealistic [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF][START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF][START_REF] Alfaro | Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates[END_REF]).

However, since we have several observations we may model the variation in this parameter with a statistical modeling approach. We explore several such approaches below. Firstly, we consider approaches based on deterministic variation between two or more groups of clades. Secondly, we consider how temporal variation in speciation and extinction rates can be modeled and assessed. Thirdly, we develop an approach handling unknown variation based on mixture modeling, including the combination of mixtures with a linear modeling of the speciation rate. Finally, we attack the problem of estimating extinction rates.

Variation Among Clades

A simple way to model variation in diversification among clades is to assume that there are two categories: some clades diversify with speciation rate λ 1 and the others with rate λ 2 . The data are made of n 1 and n 2 clades in each category, respectively. The likelihood function is:

n 1 i 1 =1 Pr(x i 1 |λ 1 ) n 2 i 2 =1
Pr(x i 2 |λ 2 ).

Note that each clade is assigned to a category a priori, although there is no assumption on whether λ 1 is greater, or smaller, than λ 2 . The null hypothesis λ 1 = λ 2 can be tested by fitting this model and the null model whose likelihood is given by equation 2: the likelihood-ratio test (LRT) comparing these two models follows a χ 2 distribution with df = 1. An alternative is to use the Akaike information criterion [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF].

The present approach is easily generalized to more than two categories: let us denote the number of categories as K, then the likelihood function would become the product of K products:

K j=1 n j i j =1 Pr(x i j |λ j ),
where n j is the number of clades in the jth category. The LRT comparing this model with the null model of homogeneous diversification follows a χ 2 with df = K -1.

These models assume, mostly for simplicity, that there is no extinction (µ = 0); however, variation in extinction rate can be incorporated in a straightforward way. For instance a model with two categories diversifying with the same λ but with different extinction rates has the following likelihood function:

n 1 Pr(x i 1 |λ, µ 1 ) n 2 Pr(x i 2 |λ, µ 2 ),
which could be compared with the null model with µ > 0 whose likelihood is:

N Pr(x i |λ, µ),
with N = n 1 + n 2 . This test is related, but not identical, to the tests of equal diversification using sister-clades where the ages of clades are not needed (Paradis 2012b).

The supplementary materials provide annotated R code explaining how to build and fit any model following the present approach.

Linear Modeling

Following the previous section, two extreme models can be defined: the simplest one where all clades diversify at the same rate, and the most complex one where each clade has its own parameter(s). This second model will be overparameterized for a likelihood approach. Nevertheless, it is possible to model variation in diversification parameters with linear models. For in-stance, we may know a priori some variables that are likely to affect the value of speciation rate (e.g., body size), and a model that relates such 'covariates' to speciation rate may be an appropriate candidate to model the variation in diversification among clades. We use here a standard strategy to model variation in a rate with respect to a covariate z:

g(λ i ) = βz i + α,
where λ i is the speciation rate in clade i, g is a function used to transform the rate in order to linearize the relationship, and β and α are two parameters.

Here β controls the effect of z on λ: if β > 0 then species with large values of z will speciate faster than those with small values of z (and inversely if β < 0). It is possible to consider more than one predictor in which case the number of parameters is equal to the number of predictors plus one.

Nonlinear models can also be considered. Each clade has its own speciation rate given by (with g -1 being the inverse transformation of g):

λ i = g -1 (βz i + α), (3)
which is used to calculate the likelihood defined by equation 2: the likelihood function is then maximized to estimate β and α (see code in the Supplementary Material). A common choice for g is the logit function, ln(λ i /(1 -λ i )), so g -1 gives:

λ i = 1 1 + e -(z i β+α) ,
The null model is defined by fixing β = 0 in which case λ = 1/(1 + e -α ) for all clades. The logit function is well suited for parameters varying between 0 and 1 which is the case for speciation rates considered on geological time scales (million of years). However, speciation rates may be larger than one on shorter scales. Other transformations can be used such as the one used below.

It must be noted that the variation among clades as modeled in the previous section is a special case of linear models where the membership of a clade to a category is coded with a discrete variable and this variable is entered as a predictor into the linear model after coding it into binary 0/1 variable(s) (see appendix in Paradis 2005, for details). Therefore, continuous and categorical predictors can be combined in the linear model. Temporal Variation [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF] studied the birth-death model in a very general way, including the cases where λ and µ vary through time. Thus it is possible to derive the probability density of the distribution of the x i 's when diversification changed through time. The likelihood can be defined and fit in the same way as above. Such a temporal model can be compared with the null model of constant diversification with a χ 2 test whose df will be equal to the number of additional parameters in the first model. As before, temporal variation may reflect speciation and/or extinction rate(s). The simplest temporal model has two rates before and after a given time point in the past, so it has one additional parameter than the null model. Note that if the time point is unknown, it could be estimated from the data so there would be two additional parameters. However, a wide variety of temporal models can be defined in ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]) using the function dbdTime where the temporal variation is defined by the user with a standard R function.

Unknown Variation

The above models assume that diversification parameters vary in relation to some known variables, either categorical or continuous. On the other hand, it is possible that these variables are not observable. Such unknown variation can be modeled with two approaches depending on whether we assume that the diversification parameters vary in a discrete or continuous manner.

A mixture of distributions is based on the assumption that observations come from two or more categories each characterized by its own distribution, but the assignment of an observation to a particular category is unknown (see [START_REF] Flury | Gender identification of water pipits (Anthus spinoletta) using mixtures of distributions[END_REF], for a biological example). As a simple example, consider a mixture of two Yule processes, then the likelihood function will be:

N i=1 f Pr(x i |λ 1 ) + (1 -f ) Pr(x i |λ 2 ), ( 4
)
where f is the proportion of clades in the first category. This model has three parameters (λ 1 , λ 2 and f ) and can be compared with the null model of homogeneous speciation with a LRT with df = 2. The idea is easily generalized to more than two mixtures: a mixture with K Yule models would have 2K-1 parameters. As above, the mixture may involve speciation and/or extinction rate(s). By contrast to the situation above where clades were assigned to categories a priori, there is here no assignment a priori.

On the other hand, assignment a posteriori is possible by calculating the relative contributions to the likelihood function.

The idea may even be further generalized to include mixtures of linear models. Suppose we know that one variable, say body size, has a significant effect on speciation rate but there is some other, unknown, variation in this parameter that we want to model with a mixture. Then it is possible to calculate the λ i 's with equation 3 and use them to compute the likelihood with eq. 4. Each category would have its own parameters β and α, so a model with K categories has 3K -1 parameters.

The second approach assumes that, in the case of a Yule model, λ varies continuously across clades following a specified distribution whose parameters are estimated from the data. A transformation of λ is useful so that it follows a normal distribution: g(λ) ∼ N (µ λ , σ 2 λ ). A useful transformation here is the complementary log-log transformation: g(λ) = ln(-ln(λ)). As above we do not know the value of λ for a given clade, but this time instead of a discrete sum we have to do a continuous integration. The likelihood function is thus:

N i=1 ∞ -∞ f N (u|µ λ , σ 2 λ ) Pr(x i |g -1 (u))du,
where f N is the density function of the normal distribution. A graphical representation of the variation in λ is obtained with the inverse transformation g -1 (u) = exp(-e u ) with the dentity of u computed with the normal distribution and the estimates μλ and σ2 λ .

Estimating Extinction Rates

The estimation of extinction rates in the absence of fossil data has appeared to be a complicated issue [START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF][START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF][START_REF] Mcpeek | The ecological dynamics of clade diversification and community assembly[END_REF][START_REF] Aldous | Five statistical questions about the tree of life[END_REF][START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF][START_REF] Didier | The reconstructed evolutionary process with the fossil record[END_REF][START_REF] Hallinan | The generalized time variable reconstructed birth-death process[END_REF]. To try to tackle this problem, we implemented a procedure which fits a birth-death model estimating λ and µ simultaneously. These estimates are denoted as λBD and μBD .

Simulation Study

The present statistical modeling approach offers many possibilities and it would take a large number of simulations to assess the statistical properties of all of them. Instead, we focus on a few key questions. What is the statistical power to detect a difference in diversification between two groups of clades? How powerful is the test to detect temporal variation in diversification? What is the statistical power to detect unknown variation in diversification between two groups of clades using mixtures? Finally, what is the precision of the λ and µ estimators?

To address these four questions, we ran four sets of simulations. First, we considered a simple two-category scenario with n 1 and n 2 clades simulated with rates λ 1 and µ 1 and λ 2 and µ 2 , respectively. The times of evolution were drawn from a uniform distribution: t i ∼ U(10, 20). A phylogeny was simulated under a birth-death process during a time t i using ape starting from a single species. The number of species surviving at time t i , x i , was extracted and the pairs (x i , t i ) were analyzed as described above using a Yule model. The LRT testing the null hypothesis of homogeneous diversification was computed, and the rejection rate was assessed under different sets of parameter values: n 1 = n 2 = {1, 3, 5, 10, 20}, λ 1 = {0.1, 0.15, 0.2}, λ 2 = 0.1,

µ 1 = {0, 0.05}, µ 2 = {0, 0.05}.
Second, we performed simulations under three scenarios with different values of diversification rates before and after 30 time units. We first generated 100 values of t from a uniform distribution between 10 and 50. We then simulated clades with constant, increasing, or decreasing diversification rate. The number of species was extracted as before, and two models were fitted: the null Yule model of constant diversification, and an alternative model assuming different speciation rates before and after 30 time units (as above µ = 0 was assumed). The rejection rates of the LRTs comparing both models were computed.

Third, a scenario similar to the first one was considered: the difference is that the simulated clades were not identified to a particular category so the data were analyzed with a mixture of Yule models. We used K = 2,

n 1 = n 2 = {10
, 20, 50}, and t i ∼ U(10, 20). Four combinations of speciation and extinction rates were used: (i) the null hypothesis is true and there is no extinction:

λ 1 = λ 2 = 0.1, µ 1 = µ 2 = 0; (ii) the null hypothesis is false
and there is no extinction:

λ 1 = 0.1, λ 2 = 0.2, µ 1 = µ 2 = 0; (iii) the null
hypothesis is false but only µ varies:

λ 1 = λ 2 = 0.2, µ 1 = 0, µ 2 = 0.1; and
(iv) same than before with stronger variation in µ:

λ 1 = λ 2 = 0.2, µ 1 = 0, µ 2 = 0.15.
Finally, we performed an assessment of the precision of the estimators of speciation and extinction rates using five combinations of λ and µ: (0.1, 0), (0.1, 0.03), (0.1, 0.06), (0.2, 0.1), and (0.2, 0.15). Here t i ∼ U(10, 30) and n = 100.

The simulations were replicated 1000 times. Annotated R (R Development Core Team 2012) code is available in the Supplementary Material with guidelines on how to run these simulations so that the readers can adapt them to their own problems. Besides, we did not attempt to compare our method with previous ones because some scenarios considered here cannot be analyzed by the latter (e.g., the third scenario does not seem to be tractable with Magallón and Sanderson's method).

Application to Fish Data

We used the data from [START_REF] Vega | Why are there so few fish in the sea?[END_REF] who compiled the number of species, stem age, and percentage of marine fish species for 22 orders and from the combined analysis of phylogeny and species richness data [START_REF] Paradis | Analysis of diversification: combining phylogenetic and taxonomic data[END_REF]. All data were unmodified from the original publication and are available at http://dx.doi.org/10.1098/rspb.2012.0075. With this data set, we explored the variation in diversification using different mixtures of Yule and birth-death models. We also tried to assess whether this variation is due to differences in the speciation or in the extinction rates.

Results

Simulation Study

The first set of simulations showed that, overall, the LRT testing for different diversification rates between two categories of clades had satisfactory statistical properties (Table 1). The type I error rate (rejection rate when the null hypothesis is true, i.e., λ 1 = λ 2 and µ 1 = µ 2 ) was, as expected, close to 5% (first and seventh lines in Table 1). However, when λ -µ was the same in both categories, the rejection rate was greater than 5% (eighth line in Table 1) showing that the present test does not test for equal diversification rate. In the cases where the null hypothesis was not true, the rejection rate varied as expected: it was greater for larger sample sizes (n 1 )

and for larger contrast in the speciation or extinction rate. Interestingly, if one category of clades had smaller µ while λ was the same, then the test was able to detect this difference; however, the statistical power was less than when the same contrast in diversification was due to different λ (compare the second and third lines in Table 1).

In the second set of simulations, the test for temporal variation rejected the null hypothesis in more than 90% when µ = 0 and λ varied, whether this was an increase or a decrease (third to sixth lines in Table 2). On the other hand, the results were contrasted when µ > 0. When there was no temporal variation in the parameters, the type I error rates were inflated in relation to the value of µ (seventh and eighth lines in Table 2). When µ varied through time, the test behaved very differently depending on the direction of this variation: it did not reject the null hypothesis in most cases when µ increased (nineth line in Table 2) while it rejected it in 68% of the cases when µ decreased (tenth line in Table 2). To further investigate this point, we repeated some of these simulations but this time the null model was a birth-death model with λ and µ constant through time, and the alternative model was with λ constant and µ allowed to vary before and after 30 time units. In this situation, the test behaved as expected: the rejection rate was less than 5% when µ was constant, whereas it varied between 8% and 31% when the null hypothesis was false (Table 3). It is noteworthy that the present test to detect time-dependent extinction rate is not very powerful: it was necessary to simulate a strong contrast in µ to reach a statistical power greater than 0.2.

The third set of simulations showed that the mixture-based LRT was able to detect heterogeneous diversification among two unknown categories of clades (Table 4). The test was more powerful when the contrast was due to different λ compared to different µ. Otherwise, the test showed satisfactory statistical performance: its power increased with sample size and/or contrast in the parameters.

The distribution of the estimates of speciation rate under the Yule model, λYule , shows that this estimator appeared unbiased when µ = 0 (Fig. 1A).

On the other hand, when µ > 0 it was negatively biased though it can be observed that λYule > λ -µ so this cannot be actually taken as an estimator of the net diversification rate. The estimator based on the birth-death model, λBD , appears less biased, even though the presence of extinctions seems to induce a slightly more dispersed distribution of the estimates (Fig. 1B).

The estimates of extinction rate based on the birth-death model, μBD , were almost unbiased (Fig. 1C).

Application to Fish Data

The fit of the Yule model to the fish data at the higher level (N = 22) resulted in a global estimate λYule = 0.058 (SE = 0.002; AIC = 456). We tried to fit a birth-death model which led to a much improved fit (AIC = 376); however, the likelihood function had a pronounced ridge on the line λ = µ (not shown). The fit of mixtures of Yule models with increasing number of categories (K) showed that the best fit was with three categories (Table 5). The parameter estimates were: λ1 = 0.041, λ2 = 0.080, λ3 = 0.013, f1 = 0.65, and f2 = 0.10. The analysis of the combined taxonomic and phylogenetic data [START_REF] Paradis | Analysis of diversification: combining phylogenetic and taxonomic data[END_REF] gave λ = 0.056 and μ = 1.83 × 10 -7 .

The analysis at the level of the families (N = 97) gave for the Yule model λYule = 0.0756 (SE = 0.0016; AIC = 1483). Like above, the fit of the birth-death model resulted in a likelihood surface with a ridge on the line λ = µ. The mixture of Yule models with the best fit had two categories (Table 5); the parameter estimates were: λ1 = 0.099, λ2 = 0.036, f = 0.42.

The analysis with a model assuming continuous variation in λ across clades gave close results for both taxonomic levels. In both cases, the model fitted well and the AIC values were smaller than for any of the previous models (Table 6). Figure 2 shows the distribution of λ inferred with the estimated parameters. Trying to introduce µ did not result in successful fits and the estimates of this parameter were close to zero. [START_REF] Vega | Why are there so few fish in the sea?[END_REF] reported the percentage of marine and freshwater species at both taxonomic levels. This was distributed very asymmetrically with most orders and families having only marine or freshwater species.

Thus, we split the data into two groups whether they had more or less than 50% of marine species. A test of different speciation rates between these groups was performed. For orders, the difference was significant (LRT:

χ 2 1 = 28.09, P < 0.001) with a larger estimate for marine orders ( λ = 0.063, SE = 0.002) compared to the freshwater ones ( λ = 0.046, SE = 0.002).

An examination of the data suggested that this result was dependent on Percomorpha which is one of the youngest clades in this data set and includes 16,625 species (Fig. 3A). Removing this clade resulted in a non-significant test (χ 2 1 = 2.10, P = 0.147, N = 21). For families, an analogous result was found with a significant test (LRT: χ 2 1 = 5.58, P = 0.018) comparing marine families ( λ = 0.079, SE = 0.002) and freshwater ones ( λ = 0.071, SE = 0.002). This result was dependent on two families older than 200 Myr (Fig. 3B): the Amiidae (one species) and Polypteridae (12). Removing these two families led to a non-significant test: χ 2 1 = 2.02, P = 0.155 (N = 95).

Discussion

The analysis of phylogenetic diversification with molecular data is enjoying a remarkable success in the literature. Some spectacular results have been accomplished using complete phylogenies (e.g., [START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF][START_REF] Hugall | Accelerated speciation in colourpolymorphic birds[END_REF][START_REF] Penney | A comparative analysis of the evolution of imperfect mimicry[END_REF]). Though complete phylogenies, possibly supplemented with fossil data, are probably the best way to investigate evolutionary diversification, the goal of our study was to show the merit of an alternative approach based on the analysis of clade data.

Our modeling approach is based on the assumption that each clade is characterized by its diversification parameters and variation among these parameters can be quantified in a statistical way. [START_REF] Bokma | Testing for equal rates of cladogenesis in diverse taxa[END_REF] and [START_REF] Paradis | Analysis of diversification: combining phylogenetic and taxonomic data[END_REF] developed a method to combine information from high-level phylogenies with clade data: both authors considered the simple constant-rate birth-death model. [START_REF] Alfaro | Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates[END_REF] used similar combined data to assess variation among clades of vertebrates using a stepwise procedure (see details in Paradis 2012a). Thus the approach in the present paper complements previous methodological developments. The possibility to quantify variation among clades with linear models seems a fruitful way to avoid overparameterization. Future applications will reinforce the relative merits of this approach.

Recently, [START_REF] Stadler | Estimating speciation and extinction rates for phylogenies of higher taxa[END_REF] developed alternative likelihood functions with respect to the way higher taxa are defined. They showed that the estimation of speciation and extinction rates vary substantially depending on these definitions. While they considered only the constant-rate birth-death model, it seems possible and interesting to include their sampling scheme into the developments presented in the present paper.

Our modeling approach ignores the background phylogeny of the clades, the set of branches that link the clades together to make a higher-level phylogeny. There are two reasons for this. First, using information from the background phylogeny is straightforward when the rates of speciation and extinction are constant and homogenous, but when this assumption is relaxed it is not simple how one must assume changes in rates in the background tree. It is clear that if a well-supported background phylogeny is available, this might give additional information which can be combined with clade data (e.g., [START_REF] Paradis | Analysis of diversification: combining phylogenetic and taxonomic data[END_REF]. However, this extra information will in most cases require its own model since it relates to older diversification events compared to clade data. On the other hand, ignoring backbone phylogeny and assuming that the clades are independent units simplifies the definition of alternative models as done in this paper. Second, though some higher-level phylogenies are available (mammals, birds), we believe these are still exceptions rather than the rule. For instance, the basal relationships of reptiles, amphibians, or fishes are still debated. Therefore, having the possibility to analyze their clade data without the need of a background phylogeny is of some general application. Furthermore, the present approach can be used when analyzing sets of clades across different phyla, for instance arthropods, echinoderms, vertabrates, etc., where the background phylogeny would not be very informative since this would branch at the origin of Metazoa.

The use of mixtures as an approach to analyze heterogeneity in diversification rates is not limited to clade data. For instance, one could model speciation and extinction rates on a fully-resolved phylogeny assuming that these parameters vary among its branches though we do not know a priori which sections of the tree evolved fast and which others evolved slowly.

Furthermore, the mixture approach can also be used to model variation in rates of trait evolution along a phylogeny. In that case, the variation may be among branches (as in the previous example), or among traits where some traits are assumed to evolve faster but we do not know which ones.

Some subtle but important facts come from the results of the simulation study. Even though most of the tests considered here assumed µ = 0, they appeared not to be tests of equal diversification. If the net diversification rates (λ -µ) were equal among clades, the tests rejected the null hypothesis in more than 5% (see eighth row of Table 1). On the other hand, if λ was equal among clades, the tests detected differences in µ. It is clear that results based only on the Yule model must be interpreted with caution.

The tests of temporal variation showed some contrasted but interesting results. When the extinction rate was zero, these tests performed very well and were able to detect either a decrease or an increase in speciation rate.

However, when extinction rate was not null, the tests based on the Yule model showed poor performance with an increased type I error rate and a high type II error rate (frequency of accepting the null hypothesis when it is false) when µ decreased through time. These poor performances were corrected if the assumtion µ = 0 was relaxed (i.e., if a null birth-death model was used in place of the Yule one), though the test had low power. Some of these results make sense: the increased type I error rate obtained with the Yule model is clearly due to the fact that a pattern of accelerated speciation can be created under a diversification process with extinction, when old lineages are mostly extinct (e.g., [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF]. On the other hand, the high type II error rate of the same model when extinction rate increased through time is somehow surprising considering the widely reported results

of slowing-down diversification (Rabosky and Lovette 2008b,a;[START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF]Etienne and Haegeman 2012, among others). Obviously, the same test was not used in these studies, so this clearly requires further investigation.

Besides, the result that the test based on a birth-death model shows statistically consistent results (i.e., the null hypothesis was rejected in less than 5% when µ was constant and in more than 5% when this parameter varied through time) is encouraging and will also be further investigated. Interestingly, this test was more powerful when the extinction rate increased through time.

A particularly interesting result comes from the precision of the estimator of extinction rate, μBD , which appears to have a very small bias, even when the data were simulated with a relatively large value of µ. This contrasts with previous studies showing that the estimator of extinction rate based on complete phylogenies is, overall, inaccurate except if it is small compared to the speciation rate [START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF][START_REF] Didier | The reconstructed evolutionary process with the fossil record[END_REF]). This result is important because several authors have cast doubt on the possibility to estimate with some precision extinction rates without fossils [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF][START_REF] Aldous | Five statistical questions about the tree of life[END_REF]).

The analysis with the fish data were essentially illustrative, but the results call for several comments. The present method seems successful in quantifying variation in diversification rates from a sample of clades. The difference in the results from both taxonomic levels makes sense since we expect more variation among families than among orders. The AIC values evidence that the model assuming continuous variation in λ across clades fits better than a model with discrete variation in this parameter. Since similar tests have not been done with other data, this clearly calls for further analyses before concluding whether diversification varies continuously or discretely across clades.

The apparent failure to estimate the extinction rate, µ, of fishes is disappointing since our simulation study showed that this parameter can be estimated correctly with the present approach. The fossil record shows many episodes of radiations, extinctions, and turn-over during the evolutionary history of fishes [START_REF] Friedman | Five hundred million years of extinction and recovery: a Phanerozoic survey of large-scale diversity patterns in fishes[END_REF]. So the reality is very different from the homogeneous scenario used in our simulations. Our results combined with previous studies (e.g., [START_REF] Aldous | Five statistical questions about the tree of life[END_REF] suggest that the estimators of µ are far more complex when rate heterogeneity is present which is likely the case with most real data set.

Vega and Wiens (2012) addressed the paradox of equivalent species di-versity between marine and freshwater fishes despite the fact that freshwater environments occupy a considerably smaller fraction of the Earth's surface than oceans. In particular they wondered whether this could be related to differences in diversification rates. Our results are in agreement with these authors' who tested their hypothesis by correlating the proportion of marine species in a clade with the method-of-moment estimator from [START_REF] Magallón | Absolute diversification rates in angiosperm clades[END_REF]. We found significant differences in λ between marine and freshwater clades from the raw data; however, the small difference in λ between both groups suggested the influence of one or two clades. Hopefully, the analysis of a more comprehensive data set with the statistical tools introduced in this paper will help to solve the paradox of less biological diversity in the ocean [START_REF] Mora | How many species are there on Earth and in the ocean?[END_REF]. λ ^Yule Percentage 0.06 0.08 0.10 0.12 0.14 λ = 0.1, µ = 0 0.06 0.08 0.10 0.12 0.14 λ = 0.1, µ = 0.03 0.06 0.08 0.10 0.12 0.14 λ = 0.1, µ = 0.06 0.06 0.08 0.10 0.12 0.14 λ = 0.2, µ = 0.1 0.06 0.08 0.10 0.12 0.14 q q q q q q q q q q q q q q q q q q q q q q 150 200 250 300

Age (Ma) Number of species A) q q

Freshwater Marine q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 50 100 150 200 250 300

Figure 1 .Figure 2 .

 12 Figure 1. Distribution of the estimates of λ and µ with (A) the Yule model ( λYule ) and (B and C) the birth-death model ( λBD and μBD) under five sets of parameters (values are given in the strips). Note the different scales of the x-axes. The vertical dotted lines indicate the values of λ (A and B) or µ (C) used in the simulation (not visible if outside the range of the x-axis). In all cases n = 100 clades.

Figure 3 .

 3 Figure 3. Number of species with respect to stem clade age for (A) orders and some superorders and (B) families of fish (the legend is the same for both plots).

Table 1 .

 1 Rejection rate for the test of equality of diversification rate between two categories with n1 and n2 (= n1) clades.

	λ 1	µ 1	λ 2 µ 2			n 1
				1	3	5	10	20
	0.1	0	0.1 0	0.044 0.060 0.066 0.054 0.056
	0.1	0	0.1 0.05 0.038 0.085 0.108 0.140 0.203
	0.15 0	0.1 0	0.094 0.177 0.232 0.418 0.711
	0.15 0	0.1 0.05 0.112 0.297 0.446 0.759 0.958
	0.2	0	0.1 0	0.174 0.497 0.707 0.943 0.998
	0.2	0	0.1 0.05 0.236 0.643 0.855 0.993 1.000
	0.1	0.05 0.1 0.05 0.040 0.054 0.048 0.054 0.061
	0.15 0.05 0.1 0	0.050 0.075 0.083 0.123 0.209
	0.15 0.05 0.1 0.05 0.069 0.129 0.201 0.387 0.653
	0.2	0.05 0.1 0	0.119 0.240 0.384 0.693 0.928
	0.2	0.05 0.1 0.05 0.143 0.407 0.618 0.878 0.995

Table 2 .

 2 Rejection rate for the test of temporal variation in diversification. The null model was a Yule model with constant rate, and the alternative model was a Yule model with λ allowed to take different values before and after 30 time units. The first two pairs of columns give the parameter values used for the simulations (Ancient and Recent: values before and after 30 time units).

Table 3 .

 3 Same than in Table2but the null model was a birth-death model with constant rates, and the alternative model was a model with λ constant and µ allowed to take different values before and after 30 time units.

		Ancient		Recent	Rejection rate
	λ	µ	λ	µ
	0.1 0.05	0.1 0.05	0.019
	0.1 0.075	0.1 0.025	0.080
	0.1 0.025	0.1 0.075	0.121
	0.1 0	0.1 0.08	0.313
	0.1 0.08	0.1 0	0.211

Table 4 .

 4 Rejection rate for the test of equality of diversification rate between two unknown categories using mixtures with n clades in each category.

	λ	µ		n
			10	20	50
	0.1	0	0.011 0.006 0.011
	(0.1, 0.2) 0	0.235 0.548 0.929
	0.2	(0, 0.1)	0.057 0.143 0.369
	0.2	(0, 0.15) 0.199 0.423 0.829

Table 5 .

 5 Results of fitting models to the fish data using mixtures of Yule processes with K from two to seven.

	K	Orders		Families
		ln L	AIC	ln L	AIC
	2 -171.207 348.414	-599.846 1205.691
	3 -158.381 326.763	-599.846 1209.691
	4 -171.207 356.414	-599.846 1213.691
	5 -171.207 360.414	-599.846 1217.691
	6 -171.207 364.414	-599.846 1221.691
	7 -171.207 368.414	-599.846 1225.691

Table 6 .

 6 Results of fitting a model of continuous variation in speciation rate across orders (N = 22) and families (N = 97) of fish.

		AIC	μλ (SE)	σλ (SE)
	Orders	320.396 1.221 (0.039) 0.163 (0.032)
	Families 1137.066 1.086 (0.026) 0.224 (0.022)

Acknowledgments

We are grateful to four anonymous reviewers, the Associate Editor, and Laure Kubatko for their constructive comments on previous versions of our manuscript. Financial support was provided by grant ANR-09-PEXT-008.


Emmanuel Paradis
email: emmanuel.paradis@ird.fr

Random phylogenies and the distribution of branching times

Keywords: diversity, extinction, LTT plot, speciation

Phylogenetic trees reconstructed without fossils have become an important source of information to study evolutionary processes. A widely used method to describe patterns of phylogenetic diversification is known as the lineages-through-time (LTT) plot. Recently, it has been shown that it is possible to predict the distribution of the branching times of a phylogeny, thus making possible to derive a theoretical LTT curve conditioned on diversification parameters. Here, I review some aspects related to this prediction showing how to derive it for any time-dependent model of diversification, as well as calculating a prediction interval around a theoretical LTT curve. The accuracy of the prediction interval was assessed with simulations using fixed or random tree sizes under constant-rate models as well as two models of time-dependent diversification. The prediction intervals were relatively narrower and more accurate for larger trees. The features of this approach are discussed as well as its potential applications.

Introduction

Biological diversity on Earth is the result of several physical, chemical, and biological forces which affect the rates of species origination and extinction. The investigation of evolutionary diversification with phylogenetic trees has received a lot of attention during the last two decades (see recent reviews in [START_REF] Pennell | An integrative view of phylogenetic comparative methods: connections to population genetics, community ecology, and paleobiology[END_REF][START_REF] Morlon | Phylogenetic approaches for studying diversification[END_REF]. In particular, a wide range of methods have been developed for the statistical analysis of phylogenies in order to estimate speciation and extinction parameters or to test hypotheses about their variation through time or among different lineages. One popular approach is a graphical method known as the lineage-through time (LTT) plot [START_REF] Nee | Tempo and mode of evolution revealed from molecular phylogenies[END_REF]. This method requires a dated phylogenetic tree (also known as chronogram) from which the dates of divergence are computed. A plot of the increase in the number of lineages (on the y-axis) with respect to time (on the x-axis) is then done to realise the LTT plot. The y-axis can be log-transformed if, for instance, the increase in the number of lineages is very strong close to present, so that the structure of the curve back in the past may be more visible. The initial motivation behind this method was to find a relationship between the shape of the curve displayed by the LTT plot and the temporal variation (or lack of) in the speciation and extinction rates of the lineage under study. Two approaches to this end can be found in the literature. The first one is based on deriving explicit equations of the distribution of the branching times of a tree with respect to the diversification parameters (e.g., [START_REF] Stadler | Lineages-through-time plots of neutral models for speciation[END_REF][START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF][START_REF] Hallinan | The generalized time variable reconstructed birth-death process[END_REF]. The second one, which has enjoyed some success in the recent literature, is to simulate many trees under a given model of diversification and then derive a null distribution of the LTT plot (e.g., [START_REF] Crisp | Explosive radiation or cryptic mass extinction? Interpreting signatures in molecular phylogenies[END_REF][START_REF] Pie | The macroevolutionary dynamics of ant diversification[END_REF][START_REF] Fordyce | Interpreting the γ statistic in phylogenetic diversification rate studies: a rate decrease does not necessarily indicate an early burst[END_REF][START_REF] Lorén | Molecular phylogenetics and temporal diversification in the genus Aeromonas based on the sequences of five housekeeping genes[END_REF].

Recently, [START_REF] Lambert | Birth-death models and coalescent point processes: the shape and probability of reconstructed phylogenies[END_REF] investigated in details general macroevolutionary models of speciation and extinction. An important result from this study is that under a homogeneous model of diversification, the distribution of branching times depends only on the rates of speciation and extinction. Therefore, the LTT plot of a phylogenetic tree completely summarises the information required to estimate these rates under the assumption of homogeneity (i.e., speciation and extinction rates are the same for all species at a given time though these parameters may vary through time) and the assumption that both rates are not correlated. On the other hand, remarkable progress has been accomplished on the simulation of phylogenetic trees under various diversification models with specified or random numbers of species [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF][START_REF] Stadler | Simulating trees with a fixed number of extant species[END_REF][START_REF] Höhna | Fast simulation of reconstructed phylogenies under global time-dependent birth-death processes[END_REF]. An issue that remains open concerns the links between these different methods and approaches, particularly whether explicit equations and simulations of phylogenetic trees lead to similar prediction of branching times distribution.

Another question is about the effect of simulating a phylogeny with a specified number of species on the distribution of branching times. Finally, the nature of the variability around the predicted distribution of branching times has not yet been assessed.

In the present paper, I review the question of predicting the distribution of branching times of a phylogeny given a model of temporal variation in speciation and extinction rates. Further elaborating on the work presented in [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF], I show how the calculations can be improved and applied to any model of time-dependent diversification. I then discuss how to derive a prediction interval around the theoretical LTT curve. I also review the methods for simulating phylogenetic trees from speciation-extinction models, and discuss the relationships among the different approaches. Using stochastic simulations, it is shown that the prediction interval of the LTT plot can be accurate in different settings of diversification parameters and tree sizes.

Methods

Speciation-extinction models and the distribution of branching times

Throughout this paper, we consider a dated (i.e., ultrametric) phylogenetic tree with N species.

We assume that this tree is rooted and fully dichotomous; consequently, it has N -1 (= m) internal nodes. The branching times, measured from the root, are denoted t 1 , . . . ,t m . By convention in this paper, the root is the origin of the time scale, so that t 1 = 0, and the present is denoted T .

Birth-death models have been widely used to analyse or simulate the diversification of clades. The basic assumption of these models is that each species is continuously exposed to speciation and extinction. The rates at which these events happen may vary through time. In such a time-dependent birth-death model, the probability that a species undergo speciation follows a function denoted λ(t), and similarly for the probability of extinction µ(t). This formulation implies no lineage-specific variation: at a given time, all species have the same probabilities of speciation and extinction (i.e., the assumption of homogeneity). In order to derive the predicted LTT curve under a time-dependent birth-death model, we use the probability that a lineage, originating from a single species at time t, is not extinct at time T [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF]:

Pr(t, n T ≥ 1) = e -ρ(t,T ) W (t) , (1)
where n T is the number of living species in the lineage at time T , and with the function ρ defined as:

ρ(t, T ) = T t µ(u) -λ(u)du,
and:

W (t) = e -ρ(t,T ) 1 + T t e ρ(t,u) µ(u)du .
We note that in eq. 1, t and T are given fixed values; only n T is a random variable.

We can now derive the distribution of branching times from the phylogeny. The probability of observing a branching event in a phylogeny, denoted as π(t), is given by the probability that two lineages originating at time t survive until present multiplied by the probability of a speciation event at t. The latter probability is given by the product of the speciation probability at t, λ(t),

with the number of species living at this time. Obviously, we do not know the latter but we can substitute it by its expectation: E(n t ) = e -ρ(0,t) . Consequently, we have:

π(t) = Pr 2 (t, n T ≥ 1)E(n t )λ(t). (2)
This density gives the relative probabilities of a branching time on the interval 0 < t < T , but this is not a probability density because this function does not integrate (since t is continuous) to one.

F (t) = t 0 π(u)du T 0 π(u)du . ( 3
)
This equation is used to derive the predicted curve of the LTT plot with given functions λ(t) and µ(t) [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF]. F (t) is an increasing function with increasing t, and we have F (0) = 0 and [START_REF] Höhna | Fast simulation of reconstructed phylogenies under global time-dependent birth-death processes[END_REF] derived an equation which is close but different from eq. 2 above.

F (T ) = 1.
However, it can be shown that both forms are equivalent (see Appendix A).

As one can see from the above formulae, the calculations involved in computing F (t) are quite heavy with several integrations which, in general, cannot be solved analytically (a trivial exception is the case with λ and µ constant). In practice, standard numerical methods for integration can be used but they are usually slow. [START_REF] Hallinan | The generalized time variable reconstructed birth-death process[END_REF] found an interesting solution to this problem: if one assumes that λ(t) and µ(t) are constant across some time intervals (i.e., they vary following some step functions), then eq. 1 can be calculated easily. This is an interesting solution because the assumption that speciation and extinction rates are constant over some short time periods (e.g., a few thousand years) may be reasonable on a macroevolutionary scale (several million years).

Here I develop some improvements on the computation of the above equations that make the calculation of F (t) faster for any function λ(t) and µ(t). First, we note that eq. 1 simplifies to:

Pr(t, n T ≥ 1) = 1 + T t e ρ(t,u) µ(u)du -1
.

We then log-transform eq. 2:

ln π(t) = -2 ln 1 + T t e ρ(t,u) µ(u)du -ρ(0,t) + ln λ(t).
Working on a logarithmic scale is less likely to inflate computing errors since we work here with sums instead of products. Note that if the antiderivates of λ(t) and of µ(t) can be found, say Λ(t)

and M(t), then ρ(t, u) can be computed directly:

ρ(t, u) = M(u) -Λ(u) -M(t) + Λ(t).
For instance, if both follow some step functions, as assumed in [START_REF] Hallinan | The generalized time variable reconstructed birth-death process[END_REF], then these antiderivates are calculated simply with rectangle areas. The two remaining difficulties are the integral of e ρ(t,u) µ(u) and the one of π(t). To compute them efficiently, the following strategy has been adopted. The function to be integrated is first evaluated on a sufficiently large number of points which are defined so that two neighbouring points are close enough so the area under the curve can be approximated with a rectangle. For instance (where δ is a small value):

x+δ x π(u)du ≈ π(x + δ/2) × δ.
This is a standard step in computing an integral with an unknown antiderivative and is called the midpoint method. The integral between 0 and t is then calculated with the sum of all the rectangle areas between these two values. We have to be careful that the error generated by the midpoint method, which may be negligible for a single small interval, will tend to accumulate when summing over many such intervals. A solution is to use Simpson's rule which is defined by:

x+δ x π(u)du ≈ π(x) + 4π(x + δ/2) + π(x + δ) 6 × δ.
Over a wider interval made of many elementary intervals, the values of all left-, mid-, and right-points of each elementary interval are summed using the same coefficients as above (1, 4, and 1). In the framework of the present application, it was found that defining 500 elementary intervals was enough, so that it is needed to calculate 999 values of π(t). Both integrals in eq. 3 are then calculated with Simpson's rule (the denominator being calculated over all elementary intervals).

The lineage-through time plot and its prediction interval

Simulations show clearly that the distribution of branching times is highly variable even with the same parameter values. So it is desirable to predict the variability around the expected LTT plot in the same way than prediction bands are derived in statistical practice [START_REF] Liu | Construction of exact simultaneous confidence bands for a simple linear regression model[END_REF]. We use here an approach introduced by Hallinan (2012) whose rationale is as follows: consider a point in time, say τ, and let us write p = F (τ). Consider now a branching time t taken at random in a tree generated with the same parameters used to compute F . Then, by the definition of a CDF, the event t ≤ τ follows a Bernoulli distribution with probability p. We can thus derive upper and lower bounds for the CDF of branching times by calculating the quantiles of the number of observed branching times for each value of t. Suppose we observe n branching times, then the number of those smaller than or equal to τ follows a binomial distribution with parameters n and p. It is therefore possible to calculate a 95% prediction interval of this distribution using the quantiles of the binomial distribution at the probabilities 0.025 and 0.975. In practice, these bounds are calculated at a number of points in the interval [0, T ] (every 0.1 time unit in this paper).

Simulating speciation-extinction trees

Simulating evolutionary data has become an important endeavour in data analysis with the constant development of Monte Carlo or Bayesian approaches. There is a long history of the simulation of trees with different purposes (see [START_REF] Paradis | Simulation of phylogenetic data[END_REF], for a brief overview). The simulation of speciation-extinction (or birth-death) trees is only a subset of these methods. Two broad categories of methods to simulate speciation-extinction trees can be defined: the homogeneous methods which are continuous in time, and the heterogeneous methods which can model any kind of variation in speciation and extinction rates (in relation to species biological traits, environment, climate, . . . ) I consider only the former category here, and distinguish three approaches to this:

1. Starting from a single species, evolution is simulated forward in time. The times of speciation and extinction events are simulated with time-to-event (or waiting time)

formulae [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF].

2. Starting from present and a given number of species, evolution is simulated backward in time. The times of speciation and extinction events are simulated with a coalescent approach [START_REF] Stadler | Simulating trees with a fixed number of extant species[END_REF].

3. Branching times are sampled randomly from the expected distribution described in the previous section. The topology of the tree is built by random successive agglomeration of species.

There can be many variants around these broad lines. For instance, the simulation can start from a (root) node instead of a single species. Another variation is to record or not the extinction events: if one is only interested in ultrametric trees (i.e., without fossils) then the extinct lineages may be simply discarded to speed-up calculations [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF].

These three approaches share in common that they simulate branch lengths in continuous time making comparisons easier because the outputs do not depend on the definition of arbitrary time steps. The first approach makes possible to simulate a tree with a fixed value of T (N is always random), while the second one can simulate trees with fixed N (T is random), and the third one can do simulations with both T and N fixed. The first approach has been used in earlier works to simulate trees with N fixed; however, [START_REF] Hartmann | Sampling trees from evolutionary models[END_REF] demonstrated that this approach results in biased simulated trees and presented a method to remedy this problem.

Various implementations of these approaches are available in different computer programs.

The first approach is implemented in ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]) which provides two algorithms (depending on whether extinct lineages are recorded or not) which simulate trees with any time-dependent variation in λ and µ (these parameters may be constant). The second approach is implemented in TreeSim [START_REF] Stadler | TreeSim: Simulating trees under the birth-death model. R package version 2.1[END_REF] which can simulate trees with constant λ and µ and the possibility to include random sampling of species as well as episodic (mass) extinctions. In order to complete this, I implemented in ape the possibility to simulate a tree with N fixed and any kind of time-dependent variation in both λ and µ (which still can be constant). This uses a slightly modified version of the algorithm given in [START_REF] Stadler | Simulating trees with a fixed number of extant species[END_REF]. The difference is in the way waiting times are generated (step 2 of the original algorithm), specifically, the formula used to generate a random time x starting from time t is here:

n t [λ(t) + µ(t)] exp -n t t x λ(u) + µ(u)du .
We keep in mind that in the framework used in this paper the time scale goes from past to present, so x < t since the algorithm goes backward in time. If λ and µ are constant, the above formula simplifies to n t (λ + µ)e n t (λ+µ)(t-x) which is the probability density function of an exponential distribution with rate n t (λ + µ), thus making possible to use efficient random number generators from common statistical software.

Finally, the third approach is implemented in TESS which can simulate trees with any time-dependent model [START_REF] Höhna | Fast simulation of reconstructed phylogenies under global time-dependent birth-death processes[END_REF].

Simulation study

A simulation study was conducted in order to assess the effect of tree size and diversification parameters on the predictions of the distribution of branching times. Two of the above simulation algorithms were used to generate phylogenetic trees. The first one is described in [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF] and simulates a tree in continuous time where λ and µ are allowed to vary with time under any user-specified model. This is a time-forward algorithm where T is fixed and N is variable (first approach of the previous section). The second algorithm is due to [START_REF] Stadler | Simulating trees with a fixed number of extant species[END_REF] and simulates trees in continuous time using a time-backward method with N fixed (second approach of the previous section). Both algorithms require as input the values of the parameters λ and µ and are implemented in the R (R Core Team, 2014) packages ape and TreeSim [START_REF] Stadler | TreeSim: Simulating trees under the birth-death model. R package version 2.1[END_REF], respectively. The third approach described in the previous section was not considered here because it samples branching times directly from the expected LTT curve, so a good agreement between the simulated trees and the predictions should be no surprise.

The predicted LTT plot and its 95% prediction interval were computed for a given set of parameter values. Ten thousand trees were simulated for each set of parameter values with both algorithms. For each tree, the LTT plot was calculated, and the agreement with the predictions were assessed in two ways. First, the proportion of points falling within the 95% prediction interval was counted. Second, the proportion of points falling above or below the predicted curve was also counted. Two sets of simulations were done. In the first set, λ and µ were kept constant with values given in Table 1. In the second set, two models of temporal variation were considered: the first one is an early burst model of speciation with decreasing λ through time and constant µ (Fig. 2a); the second one is a model of increasing extinction in recent time and constant λ (Fig. 3a).

Results

The proportion of points of the LTT plots falling within the prediction intervals depended mainly on the size of the tree (Table 1). For trees simulated with N = 10, about 69% of the points fell within the prediction intervals for all values of λ and µ. This rose to about 90% for N = 50, and 94% for N = 200. For N = 1000, this percentage was slightly larger when λ = 0.1 (94%) than when λ = 0.2 (93%). For trees generated with a random value of N, the percentage of points outside the prediction interval was related to the value of the difference λµ: when this difference was large, the percentage was above 93% whereas it was below 70% for the smallest values.

A similar pattern emerged from the proportion of points above the predicted curve (Table 2).

Except for trees simulated with N ≤ 20, the points were as likely to be above than below the predicted LTT curve. The same result was observed with random N. Figure 1 shows the predicted LTT plots under two sets of parameters: λ = 0.1, µ = 0, and λ = 0.3, µ = 0.2. The expected number of species is the same in both cases since the difference λµ is the same. As expected, the 95% prediction intervals depended on the size of the tree (shown for N = 20 and N = 200).

Figure 2a shows the modelled rates λ and µ under the scenario of early burst together with the predicted LTT curve. Figure 2b shows the LTT plot from a tree simulated with this model. By contrast to the cases with constant speciation and extinction rates, the proportion of points of the LTT plots falling within the prediction intervals varied little with respect to N and was around 90% (Table 3). The proportion of points below the predicted LTT curve was around 35% and did not vary much with tree size except when N was random where it was above 50% (Table 4).

Figure 3a shows the modelled rates λ and µ under the scenario of increasing extinction together with the predicted LTT curve. Figure 3b shows the LTT plot from a tree simulated with this model. With this scenario, the proportion of points of the LTT plots within the prediction intervals decreased slightly with increasing tree size ranging from around 92% to around 82% (Table 3). Similarly to the previous scenario, most of these points were below the predicted curve except when N was random (Table 4).

Discussion

The LTT plot is one of the most widely used methods to assess diversification from dated phylogenies. Its simplicity to run and its ability to display temporal phylogenetic patterns make it an attractive tool for evolutionists studying many groups. This method has some limitations, particularly some are intrinsic to the analysis of reconstructed phylogenies which ignore fossils (see review in [START_REF] Morlon | Phylogenetic approaches for studying diversification[END_REF]. In spite of these limitations, which are now well-known, the LTT plot has proven to be a valuable tool for investigating evolutionary diversification from molecular phylogenies (e.g., [START_REF] Moreau | Phylogeny of the ants: diversification in the age of angiosperms[END_REF][START_REF] Bininda-Emonds | The delayed rise of present-day mammals[END_REF][START_REF] Lorén | Molecular phylogenetics and temporal diversification in the genus Aeromonas based on the sequences of five housekeeping genes[END_REF].

This paper aims to review some recent developments on the analysis of branching times from phylogenetic trees with a focus on predicting their distribution. It is possible to compute a predicted LTT curve as well as a prediction interval around the predicted curve. Some improvements to compute this predicted LTT curve are presented above. The prediction interval is computed with the assumption that the branching times of the phylogeny follow a binomial distribution at any given time.

One issue for computing a predicted LTT curve is the justification of the parameter values. In the examples presented in this paper, these values were taken as fixed and given a priori, because the data were simulated so the real diversification parameter values of the phylogenies were known. With real data, the situation is different and the parameters must be estimated from the data. Two applications are possible in this context. First, it is possible to calculate an "empirical" confidence interval around the observed LTT plot with the same approach as described above: the values of p would thus be calculated from the observed cumulated proportions of number of lineages. Such a confidence interval is actually "model free" (i.e., they do not require to estimate the diversification parameters). The second application is to infer a predicted LTT curve from a fitted model and draw it together with the observed LTT plot. For instance, suppose a constant-parameter birth-death model is fitted to a tree, then the estimated values of λ and µ can be used to draw a predicted LTT curve and compare it with the actual distribution of branching times. Such an approach can help to assess the goodness of the fit of a particular model of diversification (see also [START_REF] Wollenberg | Recognizing the forest for the trees: testing temporal patterns of cladogenesis using a null model of stochastic diversification[END_REF][START_REF] Paradis | Testing for constant diversification rates using molecular phylogenies: a general approach based on statistical tests for goodness of fit[END_REF].

The match between the prediction intervals and the simulated trees showed various degrees of agreement. For the constant-parameter models, the degree of success of the prediction intervals was mainly related to the tree size (N). With N = 1000, almost 95% of the points from the LTT plots of the simulated trees were within the 95% prediction intervals. When trees were simulated with random N, the precision of the predictions were positively related to the net diversification rate (λµ), thus confirming the importance of tree size since trees simulated with a higher diversification rate are expected to have larger sizes (since the time of evolution of the simulated trees, T , was kept constant). An interesting result was that, for almost all simulations with constant parameters, the points of the LTT plots have the same probability to fall above or below the predicted LTT curve, with only a slightly higher proportion (ca. 0.55) of points below this curve for the smaller trees (see Table 2).

The success of the predictions was more contrasted when trees were simulated with the time-dependent models of diversification. Overall, the proportion of points of the LTT plots within the prediction inverval was similar to the contant-rate case with similar tree sizes.

However, this proportion did not change much with tree size and did not reach 95% like for the contant-rate case. On the other hand, the points tend to be on one side of the predicted LTT curve and slightly preferentially above it. [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF] and [START_REF] Etienne | Prolonging the past counteracts the pull of the present: protracted speciation can explain observed slowdowns in diversification[END_REF] proposed to fit models of diversification and estimate the parameters by minimising the deviation between the observed LTT curve and the one predicted from the model using a least squares criterion. It has been shown that this fitting procedure performs better than fitting the same model by maximum likelihood at least in some situations [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: a least squares approach[END_REF]. [START_REF] Hallinan | The generalized time variable reconstructed birth-death process[END_REF] questioned the use of LTT plots as a tool to assess the fit of a diversification model to a phylogenetic tree. The difficulty appears to be the autocorrelation of the individual points of an LTT plot: successive points tend to be above or below the predicted line as shown by the examples in this paper. However, the simulation results reported here suggest that, when speciation and extinction rates are constant through time, the points do not statistically tend to be either above or below the predicted line over the complete range of the plot, even within a single LTT plot. This seems to validate the use of a least squares approach to model fitting. On the other hand, when these rates varied through time, the points tended to be above the predicted LTT curve which may result in biased parameter estimation.

This will require further study.

To conlude, the distribution of branching times can be well predicted by the above formula (eq. 3) when the speciation and extinction rates are constant through time. The predictions are still accurate but less consistent for the time-dependent models considered here, which is a point requiring further study. A prediction interval around the predicted distribution (i.e., the LTT curve) can be derived and is accurate except for the smallest tree sizes. In spite of being more than two-decade old, the LTT plot is a flexible and very useful graphical tool to display the distribution of branching times, either from phylogenetic trees, or from predictions derived from speciation-extinction models.

The probability that a lineage originating at time t has a single species at time T used by Höhna is (Kendall, 1948, eq. 16):

Pr(t, n T = 1) = e -ρ(t,T ) W (t) 2 .
We now develop eq. 2 of the present paper:

π(t) = Pr(t, n T ≥ 1) 2 E(n t )λ(t) = e -ρ(t,T ) W (t) 2 e -ρ(0,t) λ(t) = e -ρ(t,T ) W (t) 2 e -ρ(t,T ) e -ρ(0,t) λ(t).
We note that:

e -ρ(t,T ) e -ρ(0,t) = e -ρ(0,T ) .

We also note that the fraction is equal to Pr(t, n T = 1), and reordering the terms in the product this gives:

π(t) = Pr(t, n T = 1)λ(t)e -ρ(0,T ) ,
which is identical to Höhna's equation and shows what is the proportionality factor. We note that the latter does not depend on t, so its value will cancel out in the ratio of integrals in eq. 3 (also in Höhna, 2013, eq. 9). The LTT plot from a phylogeny simulated with this model (black line) and the number of living lineages (grey line) from the same simulation.

Fig. 1 .

 1 Fig. 1. Lineage-through time plots (grey lines) for ten simulated trees with different values of speciation (λ) and extinction (µ) rates and tree size (N). The thick line shows the predicted LTT distribution, and the dashed lines delimit the 95% prediction intervals. The y-axes are scaled to one. F: cumulative proportions of lineages.

Fig. 2 .

 2 Fig. 2. (a) The early burst model of speciation showing the parameters (λ and µ) and the predicted LTT curve. (b) The LTT plot from a phylogeny simulated with this model (black line) and the number of living lineages (grey line) from the same simulation.

Fig. 3 .

 3 Fig. 3. (a) A model of increasing extinction in recent time showing the parameters (λ and µ) and the predicted LTT curve. (b) The LTT plot from a phylogeny simulated with this model (black line) and the number of living lineages (grey line) from the same simulation.

Figure 2 :

 2 Figure 2: (a) The early burst model of speciation showing the parameters (λ and µ) and the predicted LTT curve. (b)The LTT plot from a phylogeny simulated with this model (black line) and the number of living lineages (grey line) from the same simulation.

Table 1 :

 1 Proportion of points from the LTT plots falling within the 95% prediction interval. The simulations were replicated 10,000 times. All trees were simulated with T = 50 and n 0 = 1. N: number of species in the tree; λ: speciation rate; µ: extinction rate.

	λ	µ				N
			10	20	50	100	200	1000 random
	0.1 0	0.692 0.817 0.900 0.930 0.942 0.943 0.935
		0.05 0.691 0.818 0.900 0.927 0.939 0.945 0.847
		0.09 0.691 0.817 0.900 0.927 0.939 0.942 0.680
	0.2 0.1	0.690 0.817 0.899 0.926 0.937 0.931 0.933
		0.15 0.690 0.815 0.897 0.925 0.937 0.932 0.893
		0.19 0.691 0.818 0.899 0.925 0.938 0.932 0.787

Table 3 :

 3 Proportion of points from the LTT plots falling within the 95% prediction interval for two models of timedependent diversification. The simulations were replicated 10,000 times. T was random except when N was random where T = 50 and n 0 = 1.

	Model			N
		20	50	100	200 random
	Early burst	0.910 0.930 0.926 0.906 0.792
	Recent extinction 0.916 0.897 0.858 0.811 0.890

Table 4 :

 4 Proportion of points from the LTT plots falling below the predicted LTT line. See Table3for details.

	Model			N
		20	50	100	200 random
	Early burst	0.353 0.358 0.343 0.310 0.716
	Recent extinction 0.278 0.276 0.257 0.233 0.616

Acknowledgements

I am grateful to Tanja Stadler for clarification about some issues on tree simulation, and to two anonymous reviewers for their comments. This is publication ISEM 201x-xxx.

Appendix A Höhna (2013, eq. 8) derived the probability of speciation events in a phylogenetic tree as (using the notation of the present paper):