Keywords: diversification, extinction, key-innovation, rate-shift test, speciation

Tests of shift in diversification associated with key innovations or directional environmental change can be performed with sister-clade comparisons. This approach is attractive because it does not require detailed phylogenetic information. I propose a new likelihood ratio test based on fitting two models of diversification. I show how this test differs from a previous likelihood ratio test based on the geometric distribution. With simulations from a wide range of situations, I show that the new test performs much better than this test and the traditional test by Slowinski and Guyer. The proposed test performs at least as well as the species richness contrast test which has been proposed by several authors in four versions. A power analysis with low number of pairs of sister-clades showed that the new test could detect a shift in diversification with five or less pairs of sister-clades, while the diversity contrast test cannot detect any shift in this situation. The former appears as more powerful than the latter, and therefore is recommended when the number of pairs of sister-clades is low (less than ten). All other tests should not be used as the present study showed they lack statistical power and robustness.

The factors driving evolutionary diversification can be classified into two categories: intrinsic to the species (e.g., morphological traits, physiological features, genetic structure and variability), and extrinsic to them (e.g., climatic variables, food resources, interspecific competition). Both categories of factors clearly interact because the adaptive value of a species trait is linked to a particular set of environmental variables. In order to test the relative influence of these two categories of factors, we need to formulate predictions on the variation in diversification rates as can be observed in real data, either fossil or recent.

Generally, we may expect intrinsic factors to result in variation among clades, whereas extrinsic factors result in temporal variation. Statistical methods have been developed to test these two classes of hypotheses. The problem of testing for temporal variation in diversification has been shown recently to be particularly difficult [START_REF] Stadler | Lineages-through-time plots of neutral models for speciation[END_REF][START_REF] Rabosky | Extinction rates should not be estimated from molecular phylogenies[END_REF][START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: A least squares approach[END_REF]. The problem comes mainly from the fact that different combinations of speciation and extinction rates may result in very similar distributions of branching times, even when these two parameters are constant through time [START_REF] Paradis | Time-dependent speciation and extinction from phylogenies: A least squares approach[END_REF].

On the other hand, the statistical methods testing for inter-clade variation in diversification rate do not seem to suffer from the same problem. Some remarkable progress have been accomplished in the analysis of phylogenetic trees combined with discrete or continuous traits [START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF][START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF][START_REF] Fitzjohn | Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies[END_REF][START_REF] Fitzjohn | Quantitative traits and diversification[END_REF]. Such sophisticated methods have been used by [START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF] to investigate the evolutionary forces behind the maintenance of self-compatibility among species of Solanaceae.

In the absence of detailed fossil data, it seems very likely that methods based on phylogenies are the most informative. However, though phylogenetic and molecular data are accumulating at high pace, detailed phylogenies are currently available only for a restricted number of groups. In this perspective, statistical methods based on the comparison of species richness between sister-clades appears as a valuable approach to test hypotheses on variation in diversification rate. The basic data structure to test such hypotheses is a series of sister-clades where one clade is characterized by a trait or an environment that is absent or different in the other. The aim is to test whether this trait or environment has increased diversification rate.

In this paper, I introduce a new statistical test for the analysis of this kind of data. This test is related to the one introduced by McConway and Sims (2004); however, using simulations I show that the former is more robust and more powerful. I also compare these two tests with the classical test by [START_REF] Slowinski | Testing whether certain traits have caused amplified diversification: an improved method based on a model of random speciation and extinction[END_REF] and the species diversity contrast.

Tests of Diversification Shift

Consider the data as n pairs of sister-clades: x 1i and x 2i are the species richnesses for the ith pair, and the subscripts 1 and 2 are the two states of the variable under focus (trait, environment, or other). In the followings, and for the simplicity of the text, state 1 will be called key-innovation. Let us further denote as x i (i = 1, . . . , 2n) the species richness of any clade in the data. [START_REF] Slowinski | Testing whether certain traits have caused amplified diversification: an improved method based on a model of random speciation and extinction[END_REF] built a test based on calculating the probability of all possible data sets that are less compatible with the null hypothesis than the observed one. This probability is calculated for each pair of sister-clades. A combined probability testing the global null hypothesis (i.e., with all pairs) is calculated with Fisher's combined probabilities -2 ∑ n i ln p i : this follows a χ 2 distribution with 2n degrees of freedom. The alternative hypothesis is that diversification rate is higher in clades characterized by the key-innovation. Thus the test is a one-tailed test and it cannot detect whether diversification rate is lower in those clades. Because of its design, the Slowinski-Guyer test is particularly prone to type I error in situation of U-shaped species richnesses [START_REF] Goudet | An improved procedure for testing the effects of key innovations on rate of speciation[END_REF][START_REF] Vamosi | Endless tests: guidelines for analysing non-nested sister-group comparisons[END_REF]; see also the discussion below). [START_REF] Barraclough | Rate of rbcL gene sequence evolution and species diversification in flowering plants (angiosperms)[END_REF] developed a species diversity contrast test which can be written as:

sign(x 1i -x 2i ) ln max(x 1i , x 2i ) ln min(x 1i , x 2i ) , i = 1, . . . , n.
This contrast is computed for each pair resulting in n values on which nonparametric statistics are applied to test the null hypothesis. [START_REF] Vamosi | Endless tests: guidelines for analysing non-nested sister-group comparisons[END_REF] pointed out that this test shares its rationale and general structure with three others: the one used by [START_REF] Wiegmann | Diversification of carnivorous parasitic insects: extraordinary radiation or specialized dead end?[END_REF] who computed the logarithm of the ratio:

ln x 1i x 2i = sign(x 1i -x 2i ) ln min(x 1i , x 2i ) max(x 1i , x 2i ) ,
as well as the one by [START_REF] Barraclough | Sexual selection and taxonomic diversity in passerine birds[END_REF] with:

sign(x 1i -x 2i ) max(x 1i , x 2i ) x 1i + x 2i ,
and the test used by [START_REF] Sargent | Floral symmetry affects speciation rates in angiosperms[END_REF] with simply x 1i -x 2i . Vamosi and Vamosi (2005, p. 575) and [START_REF] Vamosi | Endless tests: guidelines for analysing non-nested sister-group comparisons. An addendum[END_REF] The [START_REF] Barraclough | Rate of rbcL gene sequence evolution and species diversification in flowering plants (angiosperms)[END_REF] version has the logarithm of the smallest value between x 1i and x 2i , thus this can lead to a division by zero if it is equal to one. In this situation, one needs to add one to all x i 's. In all the variants of the diversity contrast test, the test may be one-or two-tailed.

McConway and Sims ( 2004) devised an alternative test based on the probability distribution of species numbers as given by a model of random speciation and extinction with rates λ and µ, respectively. Let us denote this number by X, and assume for simplicity that there is no extinction (µ = 0).

Starting from a single species (X 0 = 1) the probability of X t , the number of species at time t, is given by Yule (1924, eq. 5):

Pr(X t = x|λ) = e -λt (1 -e -λt ) x-1 , (1)
This, therefore, applies to stem group richness. The Appendix shows the same equation for µ > 0.

McConway and Sims set β = 1 -e -λt in order to draw an analogy between equation (1) and the geometric distribution whose probability density function is

β(1 -β) x-1 .
∏ Pr(X t = x i |λ),
and under the alternative hypothesis:

n ∏ Pr(X t = x 1i |λ 1 ) Pr(X t = x 2i |λ 2 ).
The log-transformed likelihood functions are optimized over λ, and λ 1 and λ 2 , respectively. The LRT is computed with twice the difference in the maximum log-likelihood: it follows a χ 2 distribution with one degree of freedom. The derivation of this test assumes that diversification rates and times of origin are homogeneous across sister pairs. The robustness of the test to situations where these assumptions are not met is assessed in the next section.

Simulation Study Methods

The present simulation study aimed at contrasting the statistical properties of the four tests. The scenarios were selected in order to characterize situations of varying statistical power (probability of rejecting the null hypothesis when it is false). Following [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF], I first simulated data from a geometric distribution with parameters β 1 and β 2 (here and below, the subscript of the parameter stand for the trait state). This repeats some of the simulations done by [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF]. I only considered two scenarios: one where the null hypothesis is true (β 1 = β 2 = 0.9), and one where it is false (β 1 = 0.9, β 2 = 0.7).

The second set of simulations also repeats some from [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF]: this uses a negative binomial distribution, a situation they found to be less favorable to their test. The negative binomial distribution is a generalization of the geometric one with an additional parameter denoted as ν. I considered six scenarios with ν = 0.5, 0.9, or 1.2 combined with the two parameter sets from the previous paragraph.

The third set of simulations used a Poisson distribution with parameter γ 1 = 10 and γ 2 = 10, 15, or 20.

These three sets of simulations considered n = 35 like in [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF], and were replicated 10,000 times.

In the fourth set of simulations, species numbers were taken from phylogenies simulated from a speciation-extinction model with parameters λ 1 , λ 2 , µ 1 and µ 2 . Five series of parameters were selected: (i) the null hypothesis is true (λ 1 = λ 2 , µ 1 = µ 2 ); (ii) the null hypothesis was false by moderately increasing λ 2 ; (iii) the value of λ 2 was the same than in the previous series but µ 2 was also slightly increased; (iv) the increase in λ 2 was stronger; and (v) the same than previously with an increased value of µ 2 . All phylogenies were simulated starting 35 units of time before present. In order to test the robustness of the tests to the lack of homogeneity among lineages, two additional sets of simulations were run. In the fifth set, λ 1 was randomly drawn from a uniform distribution whose interval width was 0.1 and mean equal to the value chosen in the previous set; λ 2 was then fixed with the same contrast than previously. Another set of simulations where both λ 1 and λ 2 were random gave similar results which are reported in the Supplementary Information. Finally, in the sixth set the time of origin of each phylogeny was randomly drawn from a uniform distribution in the interval [10,50]. These last three sets were simulated with n = 10 or 20, and replicated 1000 times.

After running these six sets of simulations, the two most powerful tests were selected for a detailed estimation of their statistical power with very low n (1, . . . , 10). The data were simulated from a speciation-extinction model with the following parameters: λ 1 = 0.1, λ 2 = 0.12, 0.14, . . . , 0.2, and µ 1 = µ 2 = 0.09.

In all sets, the rejection rate was calculated as the proportion of tests where the null hypothesis was rejected at the usual rate of 5%. The rejection rate is an estimate of the type I error rate when the null hypothesis is true, and an estimate of the power of the test when the alternative hypothesis is true. All simulations and analyses were programmed in the R language (R Development Core Team, 2010) version 2.13.0, and function tests are available in ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]). An annotated R script is available in the Supplementary Information

showing how to perform these simulations, assess type I/II error rates as well as other analyses. A Sweave document (.Rnw) is also provided as an example of how to use this script to repeat the detailed power analysis with low n.

Results

The four versions of the diversity contrast test gave close results. Somehow surprisingly, [START_REF] Sargent | Floral symmetry affects speciation rates in angiosperms[END_REF] version performed the best among these four while [START_REF] Barraclough | Rate of rbcL gene sequence evolution and species diversification in flowering plants (angiosperms)[END_REF] performed the worst. The two other versions gave the same results which were intermediate. Therefore I report only the two first results here. Repeated simulation runs showed that the estimates of rejection rates were accurate to ±0.002. Thus, all proportions were rounded to the third digit. The detailed results are in the Supplementary Information.

In all scenarios, the new test and the two versions of the diversity contrast test performed much better than the two others (Table 1-4). The only situation where the new test performed poorly was when data were simulated from a negative binomial distribution with ν = 0.5: the type I error rate reached 0.2 but the McConway-Sims test performed even worst (Table 2). On the other hand, the diversity contrast test kept its type I error rate close to 0.05. However, when the null hypothesis was false with the same negative binomial distribution, the new test appeared slightly more powerful than the diversity contrast test when ν = 0.5. The power of both tests was close to one with ν = 0.9 or 1.2.

10

With data simulated from a Poisson distribution, the diversity contrast test appeared as the most powerful one, though the new test had comparable power, and close to one, for the most contrasting situation with γ 1 = 10 and γ 2 = 20 (Table 3).

With data simulated from phylogenies, only the diversity contrast tests and the new one showed good properties, the latter having higher power (Table 4).

The Slowinski-Guyer and McConway-Sims tests showed significant power only when the clades with the highest diversification had a low extinction rate compared to the speciation rate, that is λ 2 = 0.15 and µ 2 = 0.1. When µ 2 was increased to 0.12, the power of these two tests was lost even with n = 20 pairs while the power of the new test reached 0.6, and the power of the diversity contrasts test was 0.4 and 0.5 for Barraclough et al.'s and Sargent's version, respectively (Table 4).

Heterogeneity among lineages did not affect the performance of the tests (Table 5). On the other hand, heterogeneity in the origin time of the lineages affected negatively the power of all tests but the new one (Table 6). For the new test, two versions were computed here: the first one assumed that the times of origin t were unknown and set to 1000, and the second one used the simulated values (results within parentheses in Table 6). Surprisinly, the first version had greater power, though its type I error rate was slightly inflated.

Overall, the Slowinski-Guyer test appeared to perform very poorly as it failed to reject the null hypothesis in all scenarios. Thus the type I error rate was kept at a low level, but the power of the test was low as well. The only situation where some significant power was observed was simulating phylogenies with a λ 2 = 0.15 and µ = 0.1 (Table 4). In this situation, this test was more powerful than the McConway-Sims test but much less than the new test. The

McConway-Sims test showed variable performances: its power was good with data from a geometric or a negative binomial distribution; however, the type I error rate was quite high in these situations, reaching 0.38 with ν = 0.5. These results only partially replicate those from [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF] simulation study which showed no sign of inflated type I error rate in similar situations. Considering that this test appeared as not to be recommended because of its poor performance in other situations, I did not explore this issue further.

The detailed power analysis with low n therefore concerned Sargent's version of the diversity contrast test and the new test (Fig. 1). The diversity contrast test cannot detect any shift in diversification when n < 6. On the other hand, the new test can detect a shift even when n = 1: the estimated power was 0.412 for λ 2 = 0.2 with a single pair of clades, and was 0.986 with n = 5. The power of the new test was greater than 0.6 when 3 ≤ n ≤ 5 for λ 2 ≥ 0.16.

Discussion

Finding the factors promoting or preventing the diversification of evolving lineages is one of the current challenges for evolutionary biologists. At the interspecific level, this problem can be addressed with four types of data: fossil data, phylogenies of recent species, species richnesses combined with ages of stem/crown groups, and species richnesses of sister-clades. All other forms of data (e.g., raw species richnesses) have to be considered as inappropriate because they lack a temporal dimension. Among these four kinds of data, the last one must be considered the coarsest because the temporal dimension is only relative, thus limiting the inference on relative contrast in diversification. Nevertheless, the approach is interesting precisely because it does not require accurate phylogenetic information or molecular dating. Recent applications of sister-clade comparisons include a test of increased diversification among myrmechorous angiosperm clades [START_REF] Lengyel | Ants sow the seeds of global diversification in flowering plants[END_REF], and a comparison of gall-inducing and non-galling insect lineages [START_REF] Hardy | Gall-induction in insects: evolutionary dead-end or speciation driver?[END_REF]. [START_REF] Vamosi | Endless tests: guidelines for analysing non-nested sister-group comparisons[END_REF] reviewed the properties of some tests of sister-clade comparison: the Slowinski-Guyer test, the sign test, and the species diversity contrast test. They did not consider the McConway-Sims test in their study arguing that it is likely that this test shares some of the properties of the Slowinski-Guyer test. The simulation results presented in this paper confirm this suspicion since the rejection rates of these two tests varied somewhat in parallel.

A surprising result was that the difference contrast test used by [START_REF] Sargent | Floral symmetry affects speciation rates in angiosperms[END_REF] performed slightly better than the versions by [START_REF] Wiegmann | Diversification of carnivorous parasitic insects: extraordinary radiation or specialized dead end?[END_REF] and [START_REF] Barraclough | Rate of rbcL gene sequence evolution and species diversification in flowering plants (angiosperms)[END_REF]. This may come from the fact that the simple difference makes less assumption on the distribution of the diversity contrasts compared to other versions.

The simulations with heterogeneous diversification parameters showed an interesting results: all tests were robust to this heterogeneity in the sense that the rejection rates were very close to those observed with homogeneous parameters.

Overall, the new test presented here and the difference contrast test showed the most acceptable statistical properties. The advantage of the new test compared to [START_REF] Sargent | Floral symmetry affects speciation rates in angiosperms[END_REF] was a higher power when data were simulated from phylogenetic trees. On the other hand, the new test appeared to have an inflated type I error rate with data simulated from a negative binomial with a low variance (ν = 0.5); however, the type II error rate was the lowest among the tests considered here with this distribution.

One should be cautious in concluding that the new test is generally more powerful than the others because the simulated scenarios are only a small fraction of the possible ones. Therefore we cannot exclude that in some cases, not considered here, this will not be the case. However, the R script provided with this paper may help to assess this issue by simulating other specific scenarios. The detailed power analysis showed that the new test works even in situations where the diversity contrast test collapses (i.e., n ≤ 5). [START_REF] Vamosi | Endless tests: guidelines for analysing non-nested sister-group comparisons[END_REF] examined in detail a problem encountered with the Slowinski-Guyer test: if species richnesses are U-shaped (highly skewed but in a symmetric manner with respect to the trait analyzed), then this test will find a significant effect of the trait on diversification. They illustrate this fact with the following data: x 1 = {216, 64, 33, 2, 3, 1} and x 2 = {3, 1, 2, 33, 216, 64}. The Slowinski-Guyer test is significant: χ 2 12 = 22.63, P = 0.031. However, both series of species richnesses are identical though differently ordered. The explanation is that species richness differences contribute in an asymmetric way depending on their direction. For instance, the first pair (216, 3) contributes the P-value p 1 = 0.014, while the fifth pair (3, 216) contributes p 5 = 0.991. After logarithmic transformation (i.e., -2 ln p i ), these contribute 8.572 and 0.018 to the overall χ 2 . So the fifth pair contributes only to increasing the number of degrees of freedom. The McConway-Sims test is also significant: χ 2 6 = 28.277, P < 0.001. The Sargent test does not suffer from this problem because it tests for deviation from zero in the distribution of richness contrasts: it is presently non-significant either one-tailed (P = 0.542) or two-tailed (P = 1). The new test also avoids this difficulty because under the alternative hypothesis the maximum likelihood estimates will be the same ( λ1 = λ2 ). Therefore the likelihood values under both hypotheses will also be equal, leading to χ 2 1 = 0. In order to avoid the increased type I error rate related to Fisher's method of combining probabilities, the Z-test may be used instead [START_REF] Whitlock | Combining probability from independent tests: the weighted Z-method is superior to Fisher's approach[END_REF].

However, this does not solve the problem of the very low power of the Slowinski-Guyer test as was shown by a limited number of simulations (results not shown).

To conclude, previous studies have suggested to not use the Slowinski-Guyer test because of its lack of power [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF] or high type I error rate [START_REF] Vamosi | Endless tests: guidelines for analysing non-nested sister-group comparisons[END_REF]. The present study largely agrees with these conclusions. Furthermore, it is shown here that the likelihood ratio test proposed by [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF] is not robust because it lacks some statistical power when data do not follow a geometric distribution. I suggest to use either the difference contrast test used by [START_REF] Sargent | Floral symmetry affects speciation rates in angiosperms[END_REF], or the new test proposed here.

The former has an acceptable type I error rate in a wide range of situations, whereas the latter is more powerful and should be preferred when the number of pairs of sister-clades analyzed is low (n < 10) and is robust to heterogeneity in diversification rate among clades. Furthermore, the new test is the only applicable with very low n (< 6) and can have substantial power in this situation. This probability is conditioned on the lineage surviving until present (X t ≥ 1).

Setting µ = 0 leads to equation (1).

It does not appear straightforward to use these formulae to derive a test similar to the one introduced above because of the difficulty of the additional extinction parameter(s) leading to more complicated likelihood expressions.

Furthermore, there are now a more complex set of hypotheses related to shifts in diversification. The null hypothesis of no shift in diversification may be false because of an increase in speciation rate, a decrease in extinction, or both.

Table 1: Rejection rates of the null hypothesis for five tests when data were simulated from a geometric distribution with parameters β 1 and β 2 (35 pairs of sisterclades). SG: [START_REF] Slowinski | Testing whether certain traits have caused amplified diversification: an improved method based on a model of random speciation and extinction[END_REF], MS: [START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF], BHN: [START_REF] Barraclough | Rate of rbcL gene sequence evolution and species diversification in flowering plants (angiosperms)[END_REF], S: [START_REF] Sargent | Floral symmetry affects speciation rates in angiosperms[END_REF].

β 1 β 2 SG MS BHN S New test
0.9 0.9 0.0012 0.243 0.044 0.048 0.073 0.9 0.7 0.0000 0.774 0.983 0.994 1.000

Table 2: Rejection rates of the null hypothesis for five tests when data were simulated from a negative binomial distribution with parameters β 1 , β 2 , and ν. The size (or shape) parameter ν was the same in both series (35 pairs of sister-clades).

Abbreviations as in Table 1.

β 1 β 2 ν SG MS BHN S New test
0.9 0.9 0.5 0.000 0.374 0.047 0.048 0.203 0.9 0.001 0.328 0.048 0.049 0.094 1.2 0.001 0.099 0.051 0.047 0.048 0.9 0.7 0.5 0.000 0.188 0.852 0.877 0.987 0.9 0.000 0.737 0.974 0.990 0.999 1.2 0.000 0.813 0.993 0.998 1.000

  They used a likelihood function optimized over β in the case of the null hypothesis, or over β 1 and β 2 for the alternative hypothesis. The likelihood ratio is calculated for each pair of sister-clades. These ratios are then summed over the n pairs: this likelihood ratio test (LRT) follows a χ 2 distribution with n degrees of freedom. Like the diversity contrast tests, the McConway-Sims test can be one-or two-tailed.I propose to optimize the likelihood function over λ using equation (1). In theory this requires to know the dates of origin t, but in practice this appears unimportant as these can be set equal to an arbitrarily large value (say 10

3

; see simulation results below). This has an important consequence because if t is large then e -λt will tend to zero, and the transformation used by

[START_REF] Mcconway | A likelihood-based method for testing for nonstochastic variation of diversification rates in phylogenies[END_REF]

will give β ≈ 1 even if λ varies. The other major difference with the McConway-Sims test is that the 2n clades are considered as independent observations; thus a single LRT is computed. The likelihood function under the null hypothesis is: 2n

Table 3 :

 3 Rejection rates of the null hypothesis for five tests when data were simulated from a Poisson distribution with parameters γ 1 and γ 2 (35 pairs of sisterclades). Abbreviations as in Table1.

	γ 1 γ 2	SG	MS BHN	S	New test
	10 10 0.000 0.000 0.045 0.044	0.000
	10 15 0.000 0.000 1.000 1.000	0.258
	10 20 0.000 0.000 1.000 1.000	1.000

Table 4 :

 4 Rejection rates of the null hypothesis for five tests when data were simulated from a birth-death process with speciation rates λ 1 and λ 2 and extinction rates µ 1 and µ 2 . n: number of pairs of sister-clades. Abbreviations as in Table1.

	λ 1	µ 1	λ 2	µ 2	n	SG	MS BHN	S	New test
	0.1 0.09 0.1 0.09 10 0.000 0.012 0.040 0.047	0.047
					20 0.000 0.004 0.041 0.041	0.055
	0.1 0.09 0.12 0.09 10 0.015 0.023 0.118 0.126	0.190
					20 0.012 0.030 0.246 0.282	0.379
	0.1 0.09 0.12 0.1 10 0.006 0.018 0.067 0.073	0.095
					20 0.000 0.019 0.138 0.157	0.187
	0.1 0.09 0.15 0.1 10 0.259 0.129 0.443 0.515	0.722
					20 0.453 0.210 0.771 0.884	0.954
	0.1 0.09 0.15 0.12 10 0.054 0.052 0.200 0.222	0.383
					20 0.061 0.066 0.405 0.509	0.606

Table 5 :

 5 Same than in Table4except that λ 1 was randomly drawn from a uniform distribution U in the interval [0.05, 0.15] and λ 2 was fixed as λ 1 plus the value given in the table.

	λ 1	µ 1	λ 2	µ 2	n	SG	MS BHN	S	New test
	∼ U 0.09 +0 0.09 10 0.003 0.017 0.026 0.024 0.046
					20 0.000 0.009 0.044 0.052	0.044
	∼ U 0.09 +0.02 0.09 10 0.038 0.023 0.127 0.136 0.224
					20 0.049 0.033 0.233 0.297	0.382
	∼ U 0.09 +0.02 0.1 10 0.017 0.017 0.097 0.087 0.117
					20 0.013 0.024 0.125 0.150	0.179
	∼ U 0.09 +0.05 0.1 10 0.246 0.115 0.395 0.487 0.712
					20 0.419 0.198 0.768 0.864	0.947
	∼ U 0.09 +0.05 0.12 10 0.064 0.040 0.191 0.225 0.355
					20 0.121 0.048 0.403 0.513	0.635

Acknowledgments

I am grateful to Jana Vamosi and two anonymous reviewers for their constructive comments. I thank Bernard Hugueny for comments on an earlier version of the

Appendix

The probability for a lineage originated from a single species to have x species after a time t giving speciation and extinction probabilities λ and µ is [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF]:


The statistical analysis of diversification using phylogenetic trees has progressed considerably during the last decade. New insights have been given into the estimation of speciation and extinction rates [START_REF] Nee | Inferring speciation rates from phylogenies[END_REF][START_REF] Ricklefs | Estimating diversification rates from phylogenetic information[END_REF], quantifying the temporal decline in diversification following an adaptive radiation [START_REF] Rabosky | Explosive evolutionary radiations: decreasing speciation or increasing extinction through time?[END_REF], linking traits to increased or decreased diversification [START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF][START_REF] Ree | Detecting the historical signature of key innovations using stochastic models of character evolution and cladogenesis[END_REF]; [START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF][START_REF] Moore | A Bayesian approach for evaluating the impact of historical events on rates of diversification[END_REF], as well as more general issues like community assembly [START_REF] Mcpeek | The ecological dynamics of clade diversification and community assembly[END_REF]. A notable feature of this progress is the move from purely correlative methods (such as sister-group comparisons) towards approaches based on explicit models of diversification and likelihood-based inference (for a comparison, see [START_REF] Rabosky | Likelihood methods for detecting temporal shifts in diversification rates[END_REF]. In spite of these advances, two major issues remain unsolved.

First, how accurately can temporal variation in diversification be quantified? Different evolutionary mechanisms imply different temporal patterns of changes in speciation and extinction rates. For instance, an increasing speciation rate may be the result of some ecological niches being vacant or from population fragmentation, whereas an increasing extinction rate may be a consequence of high interspecific competition, habitat destruction, or decreased genetic variability. In real situations, these different mechanisms will combine to create evolutionary scenarios. For instance, an adaptive radiation may be characterized by a very high speciation rate at or near the origin of a clade that decreases subsequently [START_REF] Schluter | The Ecology of Adaptive Radiation[END_REF][START_REF] Gavrilets | Adaptive radiation: contrasting theory with data[END_REF]. On the other hand, a phylogenetic fuse process is characterized by an initial low speciation rate during some time that increases later [START_REF] Cooper | Evolutionary explosions and the phylogenetic fuse[END_REF]. A clade expands, on average, if speciation rate is greater than extinction rate, whereas it declines if this is the opposite (Fig. 1). Moreover, for a fixed difference in these rates, species compositions will show a high turn-over if both rates are high. Since the contrast between speciation and extinction rates is likely to change through time, a clade will move through different phases of expansion, decline, or high/low turn-over. These different evolutionary mechanisms may act separately or in combination, resulting in potentially complex evolutionary diversification scenarios. It appears clearly that a statistical modeling approach must take a central place in studies of diversification. Only the confrontation of alternative models will be able to test the adequacy of complex scenarios of diversification to phylogenetic data.

Second, how accurately can in diversification be disentangled in its components of speciation and extinction? This question has remained elusive since some studies emphasized the difficulties in estimating accurately extinction rate [START_REF] Kubo | Inferring the rates of branching and extinction from molecular phylogenies[END_REF][START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF]) while others emphasized the importance of considering extinction explicitly in the estimation procedure (Nee et al. 1994a;[START_REF] Nee | Inferring speciation rates from phylogenies[END_REF]. Recently, [START_REF] Quental | Diversity dynamics: molecular phylogenies need the fossil record[END_REF] demonstrated that the estimates of extinction rates in cetaceans inferred from a molecular phylogeny are greatly underestimated compared to the same estimates inferred from the fossil record. Their study is of great importance because most groups of animals and plants do not have a sufficiently documented fossil record to be able to make such a comparison between molecular and fossil data.

In one of the first contributions to the statistical analysis of diversification, Nee et al. (1994b) described how to do likelihood-based inference of temporal variation in speciation and extinction rates, but their approach does not seem to have been implemented. In a previous paper, I developed a method to test for temporal variation in diversification rate, but it can consider only a continuous increase (or decrease) through time, and speciation and extinction rates cannot be distinguished [START_REF] Paradis | Assessing temporal variations in diversification rates from phylogenies: estimation and hypothesis testing[END_REF]. [START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF] introduced a simple statistic, the γ-statistic, to test whether diversification has increased or decreased through time. This test has been widely used in the literature. Finally, [START_REF] Rabosky | Explosive evolutionary radiations: decreasing speciation or increasing extinction through time?[END_REF] implemented Nee et al.'s (1994b) likelihood approach and applied it to study decreasing diversification with several models; in particular, they considered a decline in speciation rate, an increase in extinction rate, and a combination of both. These authors suggested that this likelihood approach can be generalized to any time-dependent model, though they did not show other examples of such models. Therefore, whether we can fit more complex models of diversification is still an open question.

In this paper, I introduce several new analytical tools for the analysis of temporal variation in speciation and extinction. I derive formulae to compute the cumulative distribution function of branching times for any time-dependent birth-death model. I introduce a new method of model fitting based on least squares, and compare it with maximum likelihood. I also introduce several algorithms to simulate phylogenies under such models. This approach is illustrated with simulated and real data.

Time-Dependent Birth-Death Models

The data under consideration are a dated phylogenetic tree with N species. We assume that this tree is rooted and fully dichotomous; consequently, it has N -1 (= m) internal nodes. By convention, the root is the origin of the time scale, and the node times, measured from the root, are denoted t i (i = 1, . . . , m) with t 1 = 0, and the present is T .

Birth-death models are a class of continuous-time branching models where each species (or individual) is continuously exposed to speciation (birth) and extinction (death). The time-dependent birth-death model (or generalized birth-death model in the terminology used in [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF]) assumes that the probability of a speciation event varies with time and follows a function λ(t), and similarly for the probability of extinction µ(t). This formulation implies no lineage-specific variation: at a give time, all species have the same probabilities of speciation and extinction.

Under such a model, we can write the probability that a lineage, originating from a single species at time t, has exactly one descendant at time T :

Pr(t, n T = 1) = e -ρ(t,T ) W (t) 2 , with ρ(t, T ) = T t [µ(u) -λ(u)]du, and
W (t) = e -ρ(t,T ) 1 + T t e ρ(t,u) µ(u)du .
Similarly, the probability of a lineage being extinct during the same time interval is

Pr(t, n T = 0) = 1 - e -ρ(t,T ) W (t) ,
from which we can write the probability that a lineage is not extinct:

Pr(t, n T ≥ 1) = 1 -Pr(t, n T = 0) = e -ρ(t,T ) W (t) . (1)
For simplification, let W (t) = e -ρ(t,T ) [1 + I(t)] with: t,u) µ(u)du.

I(t) = T t e ρ(
It appears then that all the above equations depend only on two functions: ρ(t, T ) and I(t). Therefore once we have specified λ(t) and µ(t), we may derive these two functions. This is straightforward for ρ(t, T ), especially if primitives are available for λ(t) and µ(t);

otherwise this integral can be calculated numerically. For I(t), it appears that the latter method is needed in all cases. Appendix 1 shows how these functions simplify when speciation and extinction probabilities are constant through time.

In order to fit a model of time-dependent diversification, one has to specify a parametric function for λ(t) and µ(t). It is common to assume monotonically increasing or decreasing functions, in which case a logistic model is convenient, for instance for speciation:

λ(t) = 1 1 + e -(β λ t+α λ ) , (2)
where α λ and β λ are two parameters. Equation 2 is an increasing function with respect to time if β λ > 0, or decreasing if β λ < 0; if β λ = 0 it reduces to a constant probability λ = 1/(1 + e -α λ ). The function µ(t) can be written similarly with parameters α µ and β µ .

Another advantage of the logistic model is that the primitives of equation 2 can be used to accelerate the computations. However, the present approach is not limited to monotonically changing functions.

Step functions or linearly changing functions with multiple break-points are also convenient because their primitives are easily calculated with simple geometry. Furthermore these functions may lead to biologically interesting scenarios. A step model with a single break-point for speciation probability can be defined with:

λ(t) =      λ 1 t ≤ τ λ λ 2 t > τ λ .
This model has three parameters (λ 1 , λ 2 , and τ λ ). A step model with two break-points would have five parameters (see application below).

Model Fitting to the Distribution of Branching Times

Since the probabilities of observing the different elements of a phylogeny have been written down, it is possible to develop a likelihood function following the lines of Nee et al. (1994b):

ln L = 2ρ(t 1 , T ) -4 ln [1 + I(t 1 )] + (m -1)ρ(0, T ) + (m -1) ln2 + m ∑ i=2 {-2ρ(0,t i ) + ln λ(t i ) -3 ln [1 + I(t i )]} .
The details of this derivation are provided in the Supplementary Information.

It appears from this likelihood function that the function of extinction µ appears only in the integral ρ and in the product with the exponential of ρ. On the other hand, the function of speciation λ appears as a single term. I made several attempts to implement this likelihood but this always led to meaningless estimates, suggesting it is degenerate in some way or another. This observation does not prove that it is generally impossible to fit a time-dependent model by maximum likelihood, but Appendix 2 presents a theorem suggesting that it is impossible to prove that, in general, time-dependent models can be fitted by maximum likelihood.

An alternative to a likelihood approach is to consider the distribution of branching times from the phylogeny. Since at a given time t, speciation and extinction probabilities are homogeneous (they have the same values for all species), branching times summarize all the information on diversification of the clade. The probability of observing a branching event in a phylogeny, denoted by π(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 is given by the product of the speciation probability at 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 ) = 2e -ρ(0,t) . Consequently, we have:

π(t) = Pr 2 (t, n T ≥ 1)E(n t )λ(t).
We rescale these quantities with respect to t on the interval [0, T ] to obtain a cumulative density function (CDF) of the branching times:

F (t) = t 0 π(x)dx T 0 π(x)dx. (3)
It is then possible to compare this theoretical CDF with the ECDF (empirical CDF) from the tree and find the parameter values that minimize the deviation of the former on the latter using a least squares criterion to take all information along the tree into account. The ECDF of branching times is similar to a standardized lineage-through-time (LTT) plot where both xand y-scales are scaled on [0, 1]. Considering that the points of the ECDF are equally likely to be above or below the curve of the theoretical CDF, a least squares criterion is justified and makes possible hypothesis testing by comparing different models with an F-test as a ratio of sums of squares (SS) with the most parsimonious model on the numerator. The degrees of freedom of this test are given by the number of additional parameters, and the residual degrees of freedom (mtotal number of parameters). The simulation results below show that this test does not follow exactly an F distribution so that the type I error rate is inflated. So a correction is applied by dividing the ratio of sums of squares by the number of additional parameters. This has the effect of reducing statistical power, but the type I error rate is kept at a reasonable level (see below). For a given model, uncertainty of the parameter estimates can be assessed with a bootstrap, either parametric or nonparametric [START_REF] Efron | An Introduction to the Bootstrap[END_REF]. Nonparametric bootstrapping is done by sampling with replacement the observed branching times and estimating the parameters of the model for each bootstrap sample.

In the parametric version, random samples of size N are simulated from the estimated CDF. In both cases, the standard-errors of the parameters are calculated with the standard deviation of the bootstrap estimates. The proposed SS-based approach is compared to maximum likelihood below using a simulated data set.

In most studies, it is possible to consider a priori different hypotheses about the diversification of a given clade. In practice, this implies selecting different functions λ(t)

and µ(t) that quantify these different hypotheses, choosing some parameter values, calculating F (t) for these different hypotheses, and comparing them. Figure 2A shows four hypothetical scenarios of change in diversification with different patterns for speciation and extinction: some of them would result in CDF of branching times that are very close (Fig. 2B). Figure 3 shows a more striking example with two very different scenarios that would result in very similar CDFs of branching times. These results are further discussed below.

The present approach has been implemented in R (?) and the package ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]). Additional code has been developed and is available upon request.

Simulation of Phylogenies in Continuous Time

The formulae derived above allow us to write general algorithms for simulating phylogenetic trees under any time-dependent birth-death model. Three algorithms are described in Appendix C of the Supplementary Information. The general idea behind these algorithms is the following: start the simulation from the root, and consider successively whether each lineage survives until T using equation 1. If it survives, consider whether there has been a branching event on this branch before T using appropriate equations; if yes, then apply the same procedure recursively. The difference between algorithms 1 and 2 is that the latter ignores lineages that do not survive until T , thus producing a "pruned" ultrametric phylogeny (algorithm 1 produces a "complete" phylogeny). Algorithm 3 generates a random sample of branching times (i.e., without topology) from a theoretical CDF using equation 3 and the reverse method.

These algorithms are implemented in R so that the user can easily specify functions for the temporal variation in speciation and extinction rates (some examples are provided in ape). In order to assess the bias of the estimators and the type I error rate of the F-test, I simulated 1000 phylogenies with the same parameters (λ = 0.13, µ = 0.05) using algorithm 2. The number of species in these trees ranged between 9 and 879 (median = 153, mean = 181.3). Both models were fitted to each phylogeny and the F-test was computed. The observed type I error rate (proportion of tests where the null model, which is correct here, was rejected at the nominal level 0.05) was 0.083. The results are detailed in the Supplementary Information. There was no relationship between sample size (number of species) and the rate of rejection of the null hypothesis (Fig. S3). The estimator of λ showed a slight negative bias: the mean of the 1000 estimates was 0.1148 (min = 0.0417, max = 0.1937, median = 0.1157). This bias was not related to sample size; however, the variance of the estimates decreased with increasing sample size (Fig. S4). On the other hand, the estimator of µ showed systematic and considerable variance (Fig. S5). This variance decreased with increasing number of species, although the bias remained constant with a median value of μ close to zero. This appears as a discouraging result because the trees with many species did not yield more accurate estimates of µ. This is particularly contrasting with the fact λ was more accurately estimated with the large trees. The two estimates λ and μ showed strong covariation (Fig. S6). Interestingly, this line of covariation seems to converge to the point of the real parameters.

Simulation Study

Figure 5A This is reflected by wider confidence intervals for the former than for the latter (see legend of Figure 5 and details in SI). The confidence interval of αµ is remarkably large, and the estimated value is close to the upper bound of this interval. This surprising result is due to the fact that the estimation of this parameter can be greatly influenced by a few observations, so it varies accordingly among bootstrap samples. This is confirmed by the simulation results below.

As in the constant-probability case, I simulated 1000 phylogenies with the same four-parameter model and followed the protocol of model fitting and testing as above.

The number of species (N) in these simulated trees ranged between 5 and 690 (median = 115.5, mean = 138.1). The null hypothesis was rejected in 477 cases, so the power of the test was estimated to be 0.477. Surprisingly, there was no relationship between sample size and P-value: the rate of rejection of the null hypothesis was similar with small or large phylogenies (Fig. S7). The estimates of the slope of the logistic model for speciation (β λ ) showed a positive bias for the smallest trees, but this bias converged to zero with increasing N and was nearly null for N > 200 (Fig. S8). The variance of these estimates also decreased with increasing sample size. The estimates of β µ were also positively biased but with higher variance and it converged only partially to zero with increasing N; however, the variance decreased considerably for N > 150 (Fig. S9). The estimates of α µ varied greatly and their distribution was highly skewed: 996 estimated values were less than the true value (median = -7.78). This is consistent with the bootstrap results in the previous paragraph. This suggests care in the interpretation of the estimation of this parameter. Encouragingly, the bias and dispersion of these estimates seem to decrease with increasing N (Fig. S7).

Applications

Two data sets were considered: the phylogeny of Catarrhini (Old World monkeys) taken from [START_REF] Purvis | A composite estimate of primate phylogeny[END_REF], and the phylogeny of Agamidae taken from ?. These trees have 97 and 69 species, respectively. They have been chosen in order to provide relatively homogeneous clades since large phylogenies (e.g., supertrees) are likely to include several clades with different parameters of diversification [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF][START_REF] Paradis | Detecting shifts in diversification rates without fossils[END_REF][START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF]. These two clades of primates and reptiles have attracted attention from evolutionists because they are considered as adaptive radiations [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF]?).

As emphasized above, the choice of models is crucial when testing hypotheses. In the simulation study, I only considered two models because the objective was to assess the properties of the estimators and the tests. With the present applications, I fitted several models in order to quantify the temporal variations in speciation and extinction.

Simple models were fitted first, then more complicated models depending on the results.

An attempt was made to model temporal variation in extinction even though it is anticipated that this will be difficult.

Six models were fitted to the Catarrhine data set (Table 1). All tests comparing the different models with the constant-parameter model were not significant, leading to the conclusion that diversification occurred at a constant rate for these primates.

Interestingly, the model assuming that extinction probability follows a step model (4) fits better than the null model (1), though not significantly; in that case the extinction probability before 3.02 Ma was estimated to be 0.109. All models suggest a slight increase in diversification through time, though not statistically significant (Table 2). For the sake of comparison, I calculated the γ-statistic (which follows a standard normal distribution under the null hypothesis of constant diversification; [START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF] for this data set. The result was not statistically significant (γ = -0.94, P = 0.344) further indicating that the Catarrhini diversified at a constant rate.

Five models were fitted to the Agamid data set (Table 3). The two models with time-varying speciation probability fitted significantly better than the null model (P < 0.001 in both cases). Among these two models with a logistic model ( 2) or a step model ( 3), the former fitted better though it has one parameter less than the latter. None of the time-varying extinction probability models fitted better than the null model. These results provide evidence for a continuous decrease in speciation rate for the agamid lineage (Table 4). The γ-statistic also indicated that diversification rate in this lineage was not constant (γ = -4.49, P < 0.001).

Discussion

Molecular phylogenies (actually all phylogenies estimated with recent species) are valuable sources of information in the study of macroevolutionary processes. However, the fact that extinctions are not observed, and so must be inferred from the distribution of branching times is a difficulty that has remained essentially unsolved because the uncertainty caused by extinctions cannot be quantified easily [START_REF] Harvey | Phylogenies without fossils[END_REF][START_REF] Kubo | Inferring the rates of branching and extinction from molecular phylogenies[END_REF][START_REF] Paradis | Analysis of diversification: combining phylogenetic and taxonomic data[END_REF][START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF]. It is intuitive that the absence of past events implies some biases in the inference of diversification rates. Different scenarios of temporal variation in speciation and extinction will result in different expected distributions of branching times; however, some of these scenarios will be very difficult to distinguish because these distributions will be very close. One explanation for this phenomenon is that the mean values of speciation and extinction probabilities over time are relevant instead of the exact temporal variation because of the integration over time.

Nevertheless, it is not completely clear whether these different branching times distributions may be completely undistinguishable, or whether they are actually unique and could be distinguished with sufficiently large data sets. This question needs further work to be clarified.

The simulation study presented here shows some interesting results. In the case of a constant-probability model, there is a moderate bias in the estimator λ though this bias does not seem to be affected by the size of the tree. It is not clear what is the cause of this bias. The bias in the estimator μ is large and this estimator has a large variance as well. The type I error rate (probability of wrongly rejecting the null hypothesis) is low (≈ 8%). On the other hand, the same test performed on phylogenies simulated with a time-dependent model had a moderate power (≈ 50%) even though the distribution of the F-test is not completely determined. Further work on this point is clearly needed to improve the performance of the test (increase power while keeping the type I error rate at a low level). In the case of time-dependent speciation and extinction probabilities, the estimators of the slopes βλ and βµ were positively biased, but this bias converged towards zero with increasing tree size. As in the constant-probability model, the variance of the estimates related to extinction ( βµ ) were greater than for speciation. It is clear that the range of scenarios considered in this paper is restricted, and it will be necessary to assess these statistical properties for a wider range of situations.

Two previous studies showed that when speciation and extinction probabilities are both constant, the extinction probability estimator may be very imprecise [START_REF] Kubo | Inferring the rates of branching and extinction from molecular phylogenies[END_REF][START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF]). Particularly, unbiased estimation of a constant extinction probability can be achieved only when this parameter is close to zero and the age of the root of the reconstructed tree is close to the age of the actual ancestor of the lineage, in other words when the tree without dead lineages is close to the tree with them [START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF]. Recently, Rabosky (2009a[START_REF] Rabosky | Extinction rates should not be estimated from molecular phylogenies[END_REF] showed that when extinction rate is heterogeneous among species or through time some biases in the estimation of extinction rate are also observed. It must be emphasized that even when extinction rate is constant, its estimation is biased in most cases. Rabosky (2009b) developed a model of constant diversity and heritable extinction rates (further discussed below) which can produce phylogenies with estimated extinction rates equal to zero. However, this may be a special case of a more general explanation: if some species have higher extinction risk (i.e., if extinction rate is heterogeneous among species), they are less likely to be present in phylogenies of living species. Consequently, these phylogenies provide a biased sample including a majority of the species with the lowest extinction risk. This more general explanation does not require constant diversity or heritable extinction risk. Various other explanations have been proposed to explain the failure to estimate extinction rates (e.g., lack of appropriate models or of fitting methods), but it could be simply that since extinctions are not observed in phylogenies of recent species, estimating the rate at which these events occur is not generally possible.

Temporal variation in diversification parameters may create patterns of branching times similar to a constant-rate model. It is crucial to appreciate this point because many relevant hypotheses about phylogenetic diversification consider temporal variation in speciation and extinction probabilities [START_REF] Nee | Tempo and mode of evolution revealed from molecular phylogenies[END_REF][START_REF] Cooper | Evolutionary explosions and the phylogenetic fuse[END_REF][START_REF] Bininda-Emonds | The delayed rise of present-day mammals[END_REF][START_REF] Rabosky | Explosive evolutionary radiations: decreasing speciation or increasing extinction through time?[END_REF]. The striking example illustrated on Fig. 3 can be understood as both scenarios have the same mean speciation probabilities over time; consequently, the integrals computed above will be very close.

Therefore, it may be expected that complex scenarios of temporal change in diversification rates will be difficult to test. Future studies on phylogenetic diversification will need to consider this point.

There are a large number of methods to deal with model uncertainty in the statistical literature (e.g., ?Claeskens and Hjort 2008) but it does not seem easy to apply them to the present problem. Usually, techniques such as model averaging or multimodel inference are useful when the subject of interest is not the models themselves but another quantity we wish to estimate or predict. In the present context, the models are central to the question of interest since we want to know whether speciation or extinction changed through time. If we change our perspective and focus on other quantities which are products of speciation and extinction (e.g., diversity, morphological disparity, distribution, etc) then a multimodel approach might be particularly relevant because of this model indeterminacy. This raises many challenges though, such as how to tackle a potentially extremely vast model and parameter space. Some Bayesian methods may offer a way to attack this problem such as Occam's window [START_REF] Madigan | Model selection and accounting for model uncertainty in graphical models using Occam's window[END_REF] where the model space is explored in a hierarchical way.

Although maximum likelihood seems the ideal method for fitting diversification models, alternative methods have been proposed like [START_REF] Bokma | Artificial neural networks can learn to estimate extinction rates from molecular phylogenies[END_REF] neural networks. Rabosky (2009b) used approximate Bayesian computation (ABC) to fit a model of heritable extinction with pulsed turnover (HEPT). In this model, diversity is constant through time, so extinctions are balanced by speciations and the former are assumed to have correlated probabilities of occurence among species. Consequently, the HEPT model does not belong to the class of birth-death models and its likelihood function cannot be derived easily (Rabosky 2009b). However, the HEPT model is particularly interesting to model diversification because it relaxes several assumptions of birth-death models. Since the ECDF of branching times of the HEPT model can be obtained, it is possible to fit this model by the least squares method presented here and compare it to alternative models. This is a potential interesting area of future research because alternative models of diversification have already been considered in the past, but it is difficult to compare them with standard birth-death models (see, e.g., [START_REF] Losos | Stumped by trees? A generalized null model for patterns of organismal diversity[END_REF].

Incomplete taxon sampling has been of interest in the literature on diversification (e.g., [START_REF] Fitzjohn | Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies[END_REF]), particularly for extremely-rich clades where data will not be available for all species in the near future (e.g., [START_REF] Pie | The macroevolutionary dynamics of ant diversification[END_REF]. It is possible to generalize the above formulae to this situation by realizing that such subsampling of species within a clade is similar to the sampling of individuals in a population where their genealogy is modelled by the coalescent [START_REF] Kingman | On the genealogy of large populations[END_REF]. Therefore, the ECDF of branching times in the case of incomplete sampling may be obtained by superposing two processes: birth-death, forward in time, and coalescent, backward in time [START_REF] Stadler | On incomplete sampling under birth-death models and connections to the sampling-based coalescent[END_REF]. Since the number of sampled species within a clade may not be negligible compared to the total number of species, exact formulae for the coalescent should be used in this situation (see [START_REF] Wakeley | Coalescent theory: an introduction[END_REF].

The analysis of the simulated and real phylogenies showed examples with constant, declining, or increasing diversification. The results from the analysis of the Agamid data indicate that the decline in speciation probability (also found by Rabosky and Lovette 2008, with a maximum likelihood approach) was more likely to have been continuous than discontinuous as suggested by the better fit of the logistic model compared to the step model. It will be interesting to test if the same pattern can be observed in other groups since a continuous decline is more likely to be intrinsic (e.g., loss of genetic variation) whereas a discontinuous decline may be environmentally-driven.

fitting of time-dependent birth-death models can generally be done.

Model Speciation Extinction

Model Speciation Extinction

  Molecular phylogenies contribute to the study of the patterns and processes of macroevolution even though past events (fossils) are not recorded in these data. In this paper, I consider the general time-dependent birth-death model in order to fit any model of temporal variation in speciation and extinction to phylogenies. I establish formulae to compute the expected cumulative distribution function of branching times for any model, and, building on previous published works, I derive maximum likelihood estimators. Some limitations of the likelihood approach are described, and a fitting procedure based on least squares is developed that alleviates the shortcomings of maximum likelihood in the present context. Parametric and nonparametric bootstrap procedures are developed to assess uncertainty in the parameter estimates, the latter version giving narrower confidence intervals and being faster to compute. I also present several general algorithms of tree simulation in continuous time. I illustrate the application of this approach with the analysis of simulated data sets, and two published phylogenies of primates (Catarrhinae) and lizards (Agamidae). KEY WORDS: Birth-death, diversification, extinction, phylogeny, speciation.

Figure

  Figure4Ashows a tree simulated with constant probabilities λ = 0.13 and µ = 0.05. The model fitted assuming constant probabilities gave λ = 0.106 and μ = 0.004 (SS = 13.322, Fig.4B). A time-dependent model with both probabilities following equation 2 gave SS = 12.448 (Fig.4C). This second model fits better to the data but this is expected because it has two additional parameters. The test comparing them is F 2,129 = 0.535 (P = 0.587),

  Figure 5A shows a tree simulated with time-dependent probabilities according to a logistic model as shown on Fig. 5B. The model fitted with constant rates gave λ = 0.268 and μ = 0.198, and the SS was 246.654. The model with time-dependent rates gave an SS value of 15.973. The test comparing these two models is thus F 2,127 = 16.14 (P < 0.001). Figure 5B also displays the estimated curve of speciation and extinction probabilities according to the time-dependent rates model. For this particular simulated tree, I plotted the likelihood function and the SS curve around the real values of the extinction parameters (β µ = 0.02, α µ = -2). The likelihood function displayed no localor global maximum, and decreased continuously while the parameter values tend towards -∞. On the other hand, the SS curve displayed a global minimum with the parameter estimates βµ = 0.088 and αµ = -5.465 (Fig.6). The profile of the SS curve is much shallower with respect to the extinction parameters than to the speciation ones.

  109, μ2 = 0, τµ = 3.02 5 βλ = -0.014, αλ = 0.23 μ1 = 0.23, μ2 = 0, τµ = 4.02 6 λ1 = 0.268, λ2 = 0.4, λ3 = 0.358, μ = 0 τλ = 6.41, τ′ λ = 4.34

Figure 1 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 13456 Figure 1. Relationships between differences in speciation and extinction rates and phylogenetic diversification patterns.

Table 1 .

 1 Results for the Catarrhinae phylogeny. np: number of parameters; SS: sums of squares.The F -test compares the model with the constant-parameter model (1).

	Model

Table 2 .

 2 Parameter estimates for Catarrhinae. τ is the time (measured from present in Myr) when the parameter changes. The 95% confidence intervals of model 1 based on 500 bootstrap

	samples are [0.324, 0.447] and [0, 0].

Table 3 .

 3 Results for the Agamidae phylogeny. np: number of parameters; SS: sums of squares.The F -test compares the model with the constant-parameter model (1).

	Model	np	SS	F	df	P
	1 λ and µ constant	2 993.034	-
	2 λ follows a logistic model, µ constant 3	6.911 143.69 1, 67 0.0001
	3 λ follows a step model, µ constant	4	86.757	5.72 2, 66 0.005
	4 λ constant, µ follows a logistic model 3 993.000	1.00 1, 67 0.37
	5 λ constant, µ follows a step model	4 993.000	0.50 2, 66 0.61

Table 4 .

 4 Parameter estimates for Agamidae. τ is the time (measured from present in Myr) when the parameter changes. The 95% confidence intervals of model 2 based on 1000 bootstrap samples are [-0.210, 0.016], [-2.169, 0.362], and [0, 0.012].

Acknowledgments

I am grateful to Dan Rabosky for his comments on a previous version of this paper, and to two anonymous referees for their constructive comments. Financial support was provided by grant ANR-2006-BDIV002. This is publication IRD-DIVA-ISEM 2010-116.

Appendix 1

If λ and µ are constant through time we have:

From this we can derive the probability in eq. 1 when both rates are constant:

Appendix 2

Theorem. The integral I(t) has no analytical solution for any non-constant function µ(t) independent of λ(t).

Proof. To simplify the notation, functions of u are denoted with their names (i.e.,

) and all derivatives are with respect to u (i.e., ρ ′ = ∂ρ(t, u)/∂u). Suppose that I(t) has an analytical solution, therefore e ρ .µ has an antiderivative of the form f .e ρ because, from the derivative of the product of two functions, ( f .e

we should be able to find a function f such that f ′ + f .ρ ′ = µ. Since ρ ′ = µ -λ, we have:

which can be rearranged as:

Except in the case where λ and µ are constant, there will be a relationship between λ(u) and µ(u) implied by the existence of f . For instance, suppose f is a constant, so f ′ = 0 and the solution is:

This will be true only if the ratio λ(u)/µ(u) is constant over u.

This proof is only partial, but it is highly suggestive that I(t) cannot be solved analytically. Furthermore, attempts to find this integral with algebra computer programs (e.g., www.quickmath.com) failed even with some simple models such as the logistic model.

An important corollary of the above theorem is that it is impossible to prove that ML


Emmanuel Paradis
email: emmanuel.paradis@ird.fr

Jeffrey Barrett

PEGAS: an R package for population genetics with an integrated-modular approach

pegas is a new package for the analysis of population genetic data. It is written in R and is integrated with two other existing R packages (ape and adegenet). pegas provides functions for standard population genetic methods, as well as low-level functions for developing new methods. The flexible and efficient graphical capabilities of R are used for plotting haplotype networks as well as for other functionalities. pegas emphasises the need to further develop an integrated-modular approach for software dedicated to the analysis of population genetic data.

Population genetics has a strong mathematical background, and therefore genetic data analyses heavily relies on computer programs. Currently, no unified framework for these programs exists making the use of the many different population genetics programs a complicated task [START_REF] Excoffier | Computer programs for population genetics data analysis: a survival guide[END_REF]. On the other hand, R (R Development Core Team, 2009) appeared as a unified framework for analysing data in bioinformatics [START_REF] Gentleman | R Programming for Bioinformatics[END_REF] and phylogenetics [START_REF] Paradis | Analysis of Phylogenetics and Evolution with R[END_REF]. In this note, I introduce pegas (Population and Evolutionary Genetics Analysis System) which aims to fill the gap between the existing R packages and the current need in population genetic data analyses.

pegas is entirely written in R to insure maximum portability among operating systems. It requires a standard R installation as well as two packages: ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF] and adegenet [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF]. These three packages make an integrated environment for population genetic data analysis.

One of the strengths of R is the flexibility of its data structures, so that it is easy to adapt them for a particular field, resulting in much easier data input/output and manipulation (one of the most time-consuming step in data analysis). Two such structures have been derived for pegas. The main one is the class "loci" which is a simple data frame where the rows are individuals and the columns are loci and optional variables. The latter may be of any kind (continuous, integer, logical, etc). Loci are coded with R * Emmanuel.Paradis@ird.fr factors (i.e., categorical variables) which are vectors of integers. Since computing on integers is fast, counting genotypic frequencies can be done efficiently. The convention is to code a genotype with the alleles separated with a forward slash, so that once genotypes are counted, allelic frequencies are easily obtained by counting the number of alleles in each genotype. Consequently, any level of ploidy may be handled. Since most operations on allelic data require counting allele and genotype frequencies, this operation has been optimised in the function summary.loci which is used by other functions in pegas.

The class "loci" inherits the class "data.frame" (a standard class in R) so all functionalities available for the latter (manipulation, subsetting, etc.) can be used for the former [START_REF] Chambers | Software for Data Analysis: Programming with R[END_REF]. The typical example is in the use of the indexing operator [, ]. A special version of this operator has been written for the class "loci" so that the loci columns are correctly identified even after deletion of some rows and/or columns. Another advantage of the data frame structure is that each row is identified by a unique label (the rownames) so that these data may easily be matched with other data that have similar labels (e.g., a distance matrix, a DNA sequence alignment, or a genealogy).

R has many tools to analyse data frames and these may be used directly for the class "loci" and functions in pegas. A common task in population genetics is to perform some analyses for different subsets of the data (populations, breeds, plots, etc). For instance, if a data set x has a variable population, the R function by may be used to compute some statistics for each level of this variable, e.g., by(x, x$population, hw.test) to perform Hardy-Weinberg test. The second argument may be a list of variables in which case the analyses will be done for each combination of them, e.g., by(x, list(x$plot, x$treatment), summary) will compute allele and genotype frequencies for every combination of plot and treatment.

The second main special data structure of pegas is the class "haplotype" which inherits the class "DNAbin" of ape. In addition to the set of unique DNA sequences, this class includes a vector of indices identifying the individuals belonging to each haplotype. This makes possible to link a set of haplotypes with individual data (phenotypes, geographic locations, etc).

Missing data are coded explicitly in pegas. In the class "loci" the standard NA is used and these are handled in the R standard way. In the class "haplotype", the IUPAC code for ambiguous nucleotides is used and this is taken into account during calculations. The current version of pegas includes standard tools for population genetic analyses: genotypic and allelic frequencies, Hardy-Weinberg equilibrium, FST , analysis of molecular variance (AMOVA), haplotype network, mismatch distribution, Tajima's D and R2 tests for population stability, nucleotide diversity (π), the population parameter θ (= 4Neν), the site frequency spectrum, as well as several functions for reading and writing data files. pegas has also several low-level functions for manipulating its data structures (extracting observed alleled and genotypes, their frequencies, ploidy levels, building tables of all possible genotypes), making possible to extend its functionalities in a straightforward way. pegas includes basic tools for the coalescent such a computing likelihoods for a given tree, and estimating θ by maximum likelihood. All functions are accompanied with a help page describing how the calculations are done and giving the relevant literature references. These pages are compiled in a reference manual in PDF available from the above website.

There are already a large number of computer programs for handling genetic data, so it is crucial for R users to be able to analyse the data files from the most widely used ones. pegas provides a function read.loci to read allelic data in tabular form from a text file. This function has several options to specify the loci and allele sepatator as well as which columns should be considered as loci or as additional variables. pegas has several functions to convert the data structures from adegenet (class "genind") to the class "loci". Since adegenet has functions to read files created for the programs STRUCTURE, FSTAT, GENETIX, and GENEPOP (read.structure, read.fstat, read.genetix, read.genepop), these file formats may also be readily used into pegas as well. Once read into R, and possibly converted into the class "loci", data can be edited by hand with the R spreadsheet data editor (Fig. 1). pegas is distributed with a tutorial explaining, step-by-step, how to input data from different file formats (in R type vignette("ReadingFiles")).

pegas uses the efficient and flexible graphical capabilities of R [START_REF] Murrell | Chapman & Hall/CRC[END_REF]. For instance, the options for plotting a haplotype network [START_REF] Templeton | A cladistic analysis of phenotypic association with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation[END_REF] are: plot.haploNet(x, size = 1, col = "black", bg = "white", col.link = "black", lwd = 1, lty = 1, pie = NULL, labels = TRUE, scale.ratio = 1, legend = FALSE, fast = FALSE, ...)

These options make possible to plot symbols of different sizes, colours (contour and background), and the links between them may be of different colours, widths, or line types (solid, dotted, dashed, etc). These may be controlled by variables computed from the original data set with R's standard statistical and computing tools. pegas illustrates the use of an integrated-modular approach for the development of software data analysis. The three packages adegenet, ape, and pegas, complement each other for population genetics as they provide functions for spatial and multivariate analyses (adegenet), trees structures and DNA sequences manipulation (ape), and basic population genetic analyses (pegas). Furthermore, each package has enough functionalities to be used on its own, independently of the others.

ACKNOWLEDGEMENT

I am grateful to two anonymous referees for their comments on a previous version of this paper, and to Thibaut Jombart for useful discussions.

c

  Fig. 1. The R data editor showing microsatellite data. Note the first column allowing to edit the individual (row) labels.

Funding: IRD program "Spirales", and European Science Foundation (ESF) through the MolArch (Molecular Archives of Past Climatic Changes) project.

.mpl.ird.fr/pegas/

http://cran.r-project.org/web/packages/pegas/index.html