Eric Lewitus

Hélène Morlon

Detecting environment-dependent diversification from phylogenies: a simulation study and some empirical illustrations

Keywords: phylogenetics, macroevolution, diversification, environment, biodiversity

Understanding the relative influence of various abiotic and biotic variables on diversification dynamics is a major goal of macroevolutionary studies. Recently, phylogenetic approaches have been developed that make it possible to estimate the role of various environmental variables on diversification using time-calibrated species trees, paleoenvironmental data, and maximum-likelihood techniques. These approaches have been effectively employed to estimate how speciation and extinction rates vary with key abiotic variables, such as temperature and sea level, and we can anticipate that they will be increasingly used in the future. Here we compile a series of biotic and abiotic paleodatasets that can be used as explanatory variables in these models and use simulations to assess the statistical properties of the approach when applied to these paleodatasets. We demonstrate that environment-dependent models perform well in recovering environment-dependent speciation and extinction parameters, as well as in correctly identifying the simulated environmental model when speciation is environment-dependent. We explore how the strength of the environment-dependency, tree size, missing taxa, and characteristics of the paleoenvironmental curves influence the performance of the models. Finally, using these models, we infer environment-dependent diversification in three empirical phylogenies: temperature-dependence in Cetacea, δ 13 C-dependence in Ruminantia, and CO 2 -dependence in Portulacaceae. We illustrate how to evaluate the relative importance of abiotic and biotic variables in these three clades and interpret these results in light of macroevolutionary hypotheses for mammals and plants. Given the important role paleoenvironments are presumed to have played in species evolution, our statistical assessment of how environment-dependent models behave is crucial for their utility in macroevolutionary analysis.

For decades, evolutionary biologists and paleontologists have debated the relative role of major biotic and abiotic drivers of macroevolutionary dynamics (76; 84). The view that links biodiversity to abiotic change, which has been apparent to biologists since the 19th century (83; 85) and is generally referred to as the Court Jester hypothesis [START_REF] Barnosky | Distinguishing the effects of the red queen and court jester on miocene mammal evolution in the northern rocky mountains[END_REF], sees shifts in diversification as timed with climatic or geologic events. The broad-scale synchrony between global biodiversity and temperature throughout the Phanerozoic (72; 49) supports this view. More focused observations, such as the rise of angiosperms during the early Cretaceous thermal maximum (77; 18), also find support for close ties between abiotic change and species diversification. The view that links biodiversity shifts to biotic interactions, on the other hand, as outlined, for example, in Leigh Van Valen's Red Queen hypothesis [START_REF] Van Valen | A new evolutionary law[END_REF] and Ehrlich and Raven's Escape and Radiate hypothesis [START_REF] Ehrlich | Butterflies and Plants: A Study in Coevolution[END_REF], suggests that interspecific interactions (e.g., predation, mutualism, and species recognition) drive species diversification [START_REF] Liow | Red Queen: from populations to taxa and communities[END_REF]. While these various abiotic and biotic factors offer fundamentally different interpretations of how species evolve over geological time, evaluating their relative importance in practice is not trivial and has been inhibited by methodological shortcomings.

Investigations into the drivers of species diversification have been conducted using fossil specimens and phylogenetic comparative methods. Time-series of fossil specimens have been used to study speciation and extinction rates through time (70; 73), the timing and effect of mass extinction events on species richness (67; 25; 15; 46), and correlations between environmental changes and shifts in biodiversity (8; 28; 35). Fossil studies remain key to our understanding of macroevolutionary dynamics, but they are hindered by incomplete sampling [START_REF] Smith | Marine diversity through the Phanerozoic: problems and prospects[END_REF], clumping of specimens near mass extinction events [START_REF] Holland | The stratigraphy of mass extinction[END_REF], and high sensitivity to taxonomic and spatial scale (51; 26). Importantly, fossil studies are usually restricted to the few lineages for which sufficient fossil information is available [START_REF] Benton | Diversification and extinction in the history of life[END_REF].

Phylogenetic comparative methods, which tender statistical or probabilistic information from species trees, provide an alternative approach to the same type of questions, applicable to a more diverse set of organisms. They have been used to study patterns of speciation and extinction in clades (55; 52), as well as how these patterns of diversification are related to paleoenvironmental drivers (78; 87; 20; 12; 42; 16; 21; 61).

Until recently, phylogenetic approaches for analyzing the role of paleoenvironmental drivers relied on correlative techniques (20; 78; 87), often assuming that environmental changes are punctuated events. This has been improved in the last couple of years, with the development of 'environment-dependent' likelihood-based approaches [START_REF] Condamine | Macroevolutionary perspectives to environmental change[END_REF], implemented in the R package RPANDA [START_REF] Morlon | RPANDA: an R package for macroevolutionary analyses on phylogenetic trees[END_REF], that allow testing whether gradual changes in paleoenvironments had a significant influence on speciation and extinction rates, as well as quantifying the direction and magnitude of this potential influence.

Environment-dependent models take the form of time-dependent birth-death models [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF], but where the speciation and extinction rates, λ(t) and µ(t), respectively, at a given time, t, can depend on measured time-variable palaeodata, E(t). In order to test the influence of the paleoenvironment on diversification, hypotheses are made about the functional form of the dependency. Two main dependencies have been used: a linear dependency, of the form:

λ(t) = λ 0 + α λ t (1)
and

µ(t) = µ 0 + α µ t, (2)
and an exponential dependency of the form

λ(t) = λ 0 e α λ t (3)
and

µ(t) = µ 0 e αµt (4)
, where λ 0 , µ 0 , α λ , and α µ are free parameters. The probability density of observing a set of branching times in a phylogeny with a specific environmental dependence can then be computed using the likelihood expressions derived for time-dependent models, adapted to accommodate a time-dependent environmental variable (54; 19). The significance of the environmental dependence, its form (e.g. linear or exponential), its sign (positive or negative), and its strength can thus be assessed using maximum likelihood techniques.

Likelihood-based phylogenetic comparative methods have been essential aids to uncovering patterns of macroevolutionary change [START_REF] Morlon | Phylogenetic approaches for studying diversification[END_REF], and we can thus anticipate that the environment-dependent models will be of great use for understanding the relative role of major biotic and abiotic drivers on deep-time diversity dynamics. These models have already identified, for example, a positive relationship between speciation rates and global temperature in ruminants in the Cenozoic [START_REF] Cantalapiedra | Dietary innovations spurred the diversification of ruminants during the Caenozoic[END_REF], an inverse relationship between net diversification rates and temperature in bird species since the KpG boundary [START_REF] Claramunt | A new time tree reveals Earth history's imprint on the evolution of modern birds[END_REF], a positive relationship between extinction rates and sea level in birdwing butterflies [START_REF] Condamine | Deciphering the evolution of birdwing butterflies 150 years after Alfred Russel Wallace[END_REF], and a positive relationship between speciation rates and Andean uplift in neotropical orchids [START_REF] Pérez-Escobar | Recent origin and rapid speciation of neotropical orchids in the world's richest plant biodiversity hotspot[END_REF]. However, we can only have a relative confidence in these empirical results, because the performances of the environment-dependent models have not been formally assessed.

We can be confident that the likelihood expressions are correct, because the environment-dependent models are straightforward extensions of time-dependent models that have themselves been thoroughly tested using simulations (54; 31); but this does not guarantee that parameter estimation on reasonably sized phylogenies is not biased, that there is enough statistical power in the data to distinguish between alternative evolutionary scenarios (such as alternative environmental dependences), or that environmental dependence is not artifactually inferred when it did not occur. So, for example, temperature-dependency of birds diversification [START_REF] Claramunt | A new time tree reveals Earth history's imprint on the evolution of modern birds[END_REF] was inferred on a highly undersampled phylogeny (∼ 0.02%); even though the likelihoods accommodate missing taxa, can such results be trusted? Environment-dependent extinction was detected in birds [START_REF] Claramunt | A new time tree reveals Earth history's imprint on the evolution of modern birds[END_REF] and birdwing butterflies [START_REF] Condamine | Deciphering the evolution of birdwing butterflies 150 years after Alfred Russel Wallace[END_REF]. Can we trust these results given the well-known difficulty in inferring extinction rates from reconstructed phylogenies [START_REF] Morlon | Phylogenetic approaches for studying diversification[END_REF]? We cannot answer these questions, or similar questions that will undoubtedly arise in future use of environment-dependent models, unless the statistical properties of these models are thoroughly assessed using simulations.

In this paper, we compile four biotic and five abiotic paleodatasets that can potentially be used as explanatory variables in environment-dependent models and use simulated trees to evaluate the statistical performance of the corresponding models. We first focus on temperature-dependent models: we test the ability to recover speciation rates, extinction rates, and their temperature-dependencies; we also test the ability to detect temperature-dependent diversification when it occurs and check that such dependency is not wrongly detected when it does not occur. Next, because a particularly promising use of these models is to find which environmental variable, among a set, best explains diversification dynamics, we test the ability to correctly identify the environmental variable that actually influenced diversification rates and assess whether specific characteristics of the paleoenvironmental curves influence this ability. Finally, we offer best-practice guidelines for using the environment-dependent model and implement them on empirical phylogenies.

Materials and Methods

Paleodata

We collected paleodata on nine variables that have been associated with macroevolutionary hypotheses of diversification. These paleodata represent only a subset of the paleoenvironments that could potentially be used in environment-dependent diversification models, and are focused primarily on the marine environment. However, most of them are also relevant to the terrestrial environment; and they cover a wide range of biotic and abiotic hypotheses with minimal redundancy. Most importantly, they cover a wide range of resolution and temporal trends; we expect the statistical properties of the environment-dependent models on these paleodata to provide a good representation of their statistical properties in general. The paleodata can roughly be classified into three main categories. (i) Paleoclimatic variables: temperature data (inferred from δ 18 O measurements) taken from [START_REF] Zachos | An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics[END_REF]; atmospheric CO 2 data taken from [START_REF] Hansen | Climate sensitivity, sea level and atmospheric carbon dioxide[END_REF]; benthic δ 13 C data taken from [START_REF] Lazarus | Cenozoic Planktonic Marine Diatom Diversity and Correlation to Climate Change[END_REF]. Temperature is the canonical indicator of climate change. Atmospheric CO 2 varies as a result of biotic changes, such as the rise of photosynthetic plankton (60), or tectonic activity, which causes either a reduction in CO 2 caused by the subduction of carbonate-rich ocean crust or the emission of CO 2 through volcanic eruptions (5); levels of atmospheric CO 2 are thought to impact, in particular, photosynthetic organisms. δ 13 C reflects global changes in organic carbon sequestration and organic isotope fractionation ratios. In the Cenozoic, δ 13 C can represent changes in the proportion of C 4 and C 3 grasses dominating the planet, and is therefore likely to be crucial to the evolution of grasses and the herbivores that feed on them; δ 13 C can also represent increases in 13 C-rich marine species, such as diatoms, which contribute considerably to the overall carbon present in the planktonic food web (40; 59). (ii) Variables reflecting changes in oceanic composition: silica weathering ratio [START_REF] Cermeño | Continental erosion and the Cenozoic rise of marine diatoms[END_REF] and sea level estimates [START_REF] Miller | The Phanerozoic Record of Global Sea-Level Change[END_REF]. The silica weathering ratio accounts for influxes of silicic acid in the ocean, which is needed for silica-precipitating microalgae (e.g., diatoms and radiolaria), which are responsible for about a quarter of global primary production [START_REF] Buchan | Master recyclers: features and functions of bacteria associated with phytoplankton blooms[END_REF]. Sea level has several effects on both the terrestrial and the marine biome, through, for example, different overall land surfaces, different weathering levels, changes in stoichiometry, and currents. (iii) Variables reflecting interspecific interactions (primarily in the ocean): we constructed diversity curves for fossil archaeplastida, radiolaria, foraminifera, and ostracods. Fossil data were compiled from the Neptune Database [START_REF] Lazarus | Neptune: A marine micropaleontology database[END_REF] and Paleobiology Database (https://paleobiodb.org/) and diversity curves were estimated at the genus level using shareholder quorum subsampling (1) at two-million-year bins.

These diversity curves are useful for testing hypotheses of direct species interactions, such as grazing of archaeplastida by herbivores, or indirect interactions, such as competition for the same feeding resources; they may also be useful, for example, as proxies for global productivity (archaeplastida, ( 68)). Additionally, these fossil species were chosen because they are abundant enough in the fossil record to generate detailed diversity curves.

All paleoenvironmental curves were stored in the R package RPANDA freely available on CRAN [START_REF] Morlon | RPANDA: an R package for macroevolutionary analyses on phylogenetic trees[END_REF]. For the purpose of our study, to avoid biases in model selection, curves were truncated at 52 million years ago (the earliest date of the silica curve) and values were scaled to between 0 and 1.

Performance of temperature-dependent diversification models

We used simulations to test the statistical properties of the temperature-dependent diversification models. All our simulations were conducted over the time span of the environmental curves (i.e., 52 Myrs). We discarded trees with fewer than 50 tips, which always corresponded to less than 10% of the simulated trees. Throughout, speciation and extinction rates are in units of events per million years.

We tested the ability of temperature-dependent models to recover accurate parameter values. We simulated environment-dependent birth-death trees with speciation rates, λ, and/or extinction rates, µ, that varied as an exponential function of temperature T, such that at each given point in time, t, λ(t) = λ 0 e α λ T (t) or λ(t) = λ 0 and µ(t) = µ 0 e αµT (t) or µ(t) = µ 0 . For trees with a λ-dependency on T (hereafter referred to as temperature-dependent λ trees), we simulated with λ 0 = 0.1, 0.15, 0.2, fixing α λ = 0.02 and µ 0 = 0.01; α λ = -0.5, -0.2, 0.4, 0.8, 1.6 fixing λ 0 = 0.1 and µ 0 = 0.01; and µ 0 = 0.01, 0.03, 0.06, 0.09 fixing λ 0 = 0.1 and α λ = 0.02. We furthermore simulated temperature-dependent λ trees with a considerably higher µ (µ 0 = 0.7) to reflect estimates of µ inferred from the fossil record [START_REF] Foote | Principles of Paleontology[END_REF]. To do so, we increased λ too (λ 0 = 0.7, α λ = 0.35), otherwise it was difficult to simulate surviving trees. For trees with a µ-dependency on T (hereafter referred to as temperature-dependent µ trees), we simulated with λ 0 = 0.2, 0.25, 0.3, fixing µ 0 = 0.02 and α µ = 0.02; µ 0 = 0.02, 0.06, 0.1, fixing λ 0 = 0. or µ on T ) by maximum-likelihood using the RPANDA function fit_env (53) and recovered parameter estimates for λ 0 , α λ , µ 0 , and α µ . The likelihood implemented in this function is the one described in ( 54) and ( 19); here we conditioned the likelihood on stem age. Throughout, we used a Nelder-Mead optimization algorithm [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] to fit models with a dependence of λ on time or temperature and a simulated-annealing optimization algorithm [START_REF] Belisle | Convergence Theorems for a Class of Simulated Annealing Algorithms on d[END_REF] to fit models with a dependence of µ on time or temperature. This is because for α > 0.1 in the latter models the Nelder-Mead algorithm did not converge -specifically, the optimization function could not be computed at the initial parameters for a wide range of initial parameters.

Because many empirical phylogenies are undersampled, we tested the effect of undersampling on inferring parameters in temperature-dependent trees. We simulated 1000 temperature-dependent λ trees (λ(t) = 0.2e 0.4•T (t) ) and constant extinction (µ = 0.05) and 1000 temperature-dependent µ trees (µ(t) = 0.05e 0.1•T (t) ) and constant speciation (λ = 0.5).

We jackknifed the trees (i.e., sampled without replacement) at 10% intervals from 90 -50% and inferred parameter estimates by fitting a temperature-dependent speciation model to each set of 1000 trees using the function fit_env,while accounting for undersampling by specifying the corresponding sampling fraction (54; 19). The median tree size (of trees with more than 50 tips) of the original temperature-dependent λ trees was 1209 tips and of the original temperature-dependent µ trees was 1119 tips.

We tested the ability to detect temperature-dependent diversification when it occurs and to not detect it when it does not occur. Many paleoenvironmental curves, despite fluctuations in their trends, sport a general tendency to increase or decrease over time.

Therefore, in addition to testing the ability to distinguish between temperature-dependent and null, constant rate models, we tested the ability to distinguish between temperature-dependent and time-dependent models. We simulated birth-death trees with a positive and negative exponential dependency of λ on temperature (λ = 0.2e ±0.9•T , µ(t) = 0.05), with an increasing and decreasing dependency of λ on time (λ = 0.2e ±0.2•t , µ(t) = 0.05), and with constant rates (λ = 0.2, µ = 0.02). Additionally, because unaccounted-for background variation in diversification rates across lineages can potentially lead to biased inferences (64; 54; 66), we tested whether rates that vary among lineages can be misleadingly interpreted as environmental dependence. We therefore simulated trees with a 10% probability of a rate-shift in λ at each branching event, with an initial λ 0 = 0.2;

at each event where a shift occurred, the new λ was drawn from a normal distribution of values (with a standard deviation of 0.1) centered on the ancestral value; µ was held constant at 0.02. For each scenario we simulated 5000 trees; the median tree size (of trees with more than 50 tips) was 455 tips. Finally, to test whether high relative extinction rates affect our ability to properly infer temperature-dependence on speciation rates, we simulated 1000 temperature-dependent λ trees (median tree size of 190 tips) with λ = 0.75e 0.35•T (t) and µ = 0.7. Time-dependent trees were simulated using tess.sim.age [START_REF] Höhna | Fast simulation of reconstructed phylogenies under global time-dependent birth-death processes[END_REF] and trees with rate-shifts were simulated with our own code. We fit all three models (exponential dependency on temperature, exponential dependency on time, and constant rates) to each tree and compared support using the corrected Akaike Information Criterion (AICc) [START_REF] Burnham | Model Selection and Multimodel Inference[END_REF]. Fits of constant and time-dependent models were performed using the fit_bd function in RPANDA. We then calculated the percentage of each set of trees best fit by each model (the model with lowest AICc was selected). We likewise simulated birth-death trees with a fixed λ and a positive and negative exponential dependency of µ on temperature (λ 0 = 0.4, µ 0 = 0.1, α µ = ±0.4), an increasing and decreasing dependency of µ on time (λ 0 = 0.5, µ 0 = 0.01, α µ = ±0.4), and with no dependency of µ on time (i.e., constant λ 0 = 0.2, µ 0 = 0.01). For each scenario we simulated 5000 trees; the median tree size (of trees with more than 50 tips) was 358 tips. Here, in addition to fitting the corresponding models, we also fit models with an exponential dependency of λ on temperature. This was done to see whether a µ-dependency on temperature left a detectable signal in λ.

In order to have an idea of the effect of tree size on parameter estimation and the ability to recover temperature dependency, we simulated 5000 temperature-dependent λ trees (λ = 0.2e ±0.05•T (t) and µ = 0.05). We fitted the corresponding temperature-dependent model and recovered parameter estimates for each tree as above. We also fitted the exponential time-dependent and constant rate models and computed the Akaike weights of these models on each tree. Finally, we reported the results by tree size bin. This procedure can potentially introduce a bias, because trees falling in a particular size-bin are trees that, by chance, even though they share the same simulated parameter values, diversified more (or less) than others. We did not observe such a bias, however. We similarly tested the effect of tree size on parameter estimation and model recovery for temperature-dependent µ trees (µ = 0.2e 0.1•T (t) ); in these simulations, we varied λ 0 across simulations so as to avoid the bias mentioned above. λ 0 values were 0.1, 0.2, 0.3, 0.4, and 0.5.

There will necessarily be a loss of resolution between the actual fluctuations in an environmental variable over time and the reconstruction of those fluctuations with sampled data. It is therefore important to understand how sensitive the environment-dependent model is to the reconstructed fluctuations. To evaluate this, we tested the effect of discrepancies in resolution between the temperature curve used to simulate trees and the one used to fit the temperature-dependent model. We simulated pure-birth (µ = 0) trees with an exponential dependency between λ and temperature, where the smoothing function for the temperature curve was determined using generalized cross-validation [START_REF] Golub | Generalised cross-validation as a method for choosing a good ridge parameter[END_REF].

We simulated 5000 trees with positive temperature-dependency (λ = 0.2e 0.05•T (t) ) and 5000 trees with negative temperature-dependency (λ = 0.2e -0.05•T (t) ). The median tree size (of trees with more than 50 tips) was 610 tips. We then fitted temperature-dependent models to the simulated trees, using increasingly smoothed temperature curves. Smoothed curves are meant to represent degraded environmental data and were obtained by reducing the number of degrees of freedom (DOFs) applied to the smoothing function. For each tree, we compared parameter estimates, as well as statistical support in comparison with time-dependent models (with exponentially trending λ) obtained when fitting temperature-dependent models with increasingly smoothed temperature curves.

Performance of model selection on trees simulated with various environmental dependencies

Generalizing the analyses performed for temperature-dependency, we assessed whether environmental dependence can be detected on trees simulated with λ dependencies on specific paleoenvironment curves. Because we did not see any effect of extinction on the ability to recover the environment-dependent model in the case of temperature-dependency, we simulated pure-birth trees (µ = 0) with exponential dependencies between λ and each paleoenvironmental curve X (see Paleodata above). We simulated trees at various values of α λ (0.15, 0.3,0.45, 0.6, 0.75, 0.9, 1.05, 1.

to test how the strength of the simulated dependency influenced the recoverability of the environmental model. For each value of α λ and for each paleoenvironmental curve we simulated 1000 trees with a median tree size (for trees with more than 50 tips) of 310 tips.

λ 0 was identical across paleoenvironmental curves, but varied between α λ values, in order to maximize the number of surviving, reasonably sized simulated trees. For each tree, we fitted the same environment-dependent model that was used to simulate that particular tree, as well as a constant-rate λ model and a time-dependent λ model, and compared their statistical support with AICc. We selected the best model as the one with the lowest AICc support. We also compared the relative support of the models using Akaike weights.

We carried out the same type of analyses to assess whether an environmental dependence on µ can properly be inferred. We simulated trees with an exponential dependency between µ and a paleoenvironmental curve X, holding λ constant, for different values of α µ as above. For each value of α µ and for each paleoenvironmental curve we simulated 1000 trees; the median tree size (of trees with more than 50 tips) was 320. For each tree, we fitted the same environment-dependent model that was used to simulate that particular tree (µ = µ 0 e αµX(t) , λ = λ 0 ) as well as a constant-rate birth-death model and a model with µ exponentially dependent on time and constant λ. We also fitted a model where the environment influences speciation rather than extinction rates (λ = λ 0 e α λ X(t) , µ = µ 0 ). We did this under the hypothesis that trees simulated under a µ dependency on a paleoenvironmental curve could manifest a signature of that dependency in their λ.

Environmental curves differ widely in their shape and characteristics (Fig 3B). In particular, the biotic curves that we considered here (constructed from fossil data) appeared coarser than the abiotic ones. To assess whether environmental dependency was more easily detected for certain types of curves than others, we characterized the curves based on defining attributes, such as biotic/abiotic, linearity/non linearity, and six further metrics, two capturing temporal autocorrelation, two capturing global trend, and two capturing the composition of the data. Next, we assessed the effect, if any, these features had on the rate of recovery of the models at different values of α λ . To test for linearity, we tested for the presence of one or two breakpoints in the data using the threshold test [START_REF] Hansen | Testing for Linearity[END_REF].

We used two metrics reflecting temporal autocorrelation; first we measured the correlation between values in the original time-series and a new series that lags behind by some amount of time, for a series of time lags. We then fit regressions between the various time lags and the corresponding correlation value. We used the intercept of the regression to measure short-term autocorrelation, and the slope to measure the rate at which autocorrelation decreases. A global trend analysis was conducted by fitting a linear model using generalized least squares to each curve and computing the correlation, slope, and degrees of freedom for the fitted model. Finally, we estimated the average rate of change of the curves by calculating the mean over all time periods of the slope between values within each time period.

Finally, for trees simulated under each paleoenvironmental curve with an exponential dependency on λ, we compared fits for models using all the paleoenvironmental curves. For each set of simulated trees, we calculated the percentage of trees best fit (i.e., lowest AICc) by each paleoenvironmental model and reported this percentage as a function of α.

Testing macroevolutionary hypotheses of environment-dependency on empirical phylogenies

We used time-calibrated molecular trees for Cetacea, Ruminantia, and Portulacaceae. The Cetacean phylogeny was 98% sampled and constructed from six mitochondrial and nine nuclear genes in a Bayesian framework; it was time-calibrated using fossil data; and divergence dates were estimated using a relaxed molecular clock [START_REF] Steeman | Radiation of Extant Cetaceans Driven by Restructuring of the Oceans[END_REF]. For Ruminantia we used two different trees: one was constructed using full mitochondrial genomes and 16

fossil calibration points in a Bayesian framework, but was only 65% sampled (9); the other one was constructed using a supermatrix of 124 previously published trees and was 100% sampled [START_REF] Cantalapiedra | Dietary innovations spurred the diversification of ruminants during the Caenozoic[END_REF]. The two phylogenies differed markedly in their datation; in particular, the phylogeny of ( 9) hypothesized significantly more recent crown ages for many families and for Ruminantia than that of [START_REF] Cantalapiedra | Dietary innovations spurred the diversification of ruminants during the Caenozoic[END_REF]. For Portulacaceae, the tree was 40% sampled and constructed separately using a combined matrix of genomic markers for chloroplast in a maximum likelihood framework and Bayesian inference; divergence times were estimated using a relaxed clock calibrated on relevant geological events. [START_REF] Ocampo | Molecular phylogenetics of suborder Cactineae (Caryophyllales), including insights into photosynthetic diversification and historical biogeography[END_REF]. We conducted analyses on Bayesian posterior trees for Cetacea and both Ruminantia trees; we analysed both the maximum-likelihood and Bayesian tree for Portulacaceae. The above trees were chosen because diversification in each has previously been associated with an abiotic process (78; 58; 12).

The nine paleoenvironment models plus a constant-rate birth-death model and time-dependent model (with and without extinction) were fit to the empirical phylogenies, accounting for missing taxa by applying the relevant sampling fraction (54; 19). For each paleoenvironment model, we tested λ and µ as constant, linear, and exponential functions of the paleoenvironmental curves, as well as µ = 0. The best-fit model for each phylogeny was determined by ∆AICc as above, and the corresponding speciation and extinction parameter estimates were recorded. All palaeoenvironmental curves were tested. While some variables (temperature and sea level) are directly biologically relevant for the three clades, other variables (CO 2 , δ 13 C, silica, archaeplastida diversity, ostracod diversity) are directly relevant for some clades but not all, and others still (foraminifera and radiolarian diversity) are not directly biologically relevant to any of the clades. When environmental variables were not directly relevant, we still included them as negative controls. We also note that while significant support for environment-dependent models is typically interpreted by an effect of the environment on clade dynamics, such support could also occur if the rise or demise of a given clade impacts the targeted environmental variable (e.g. silica levels in the ocean are often attributed to the expansion of land plants). Hence, environmental variables that are biologically relevant for a given clade are those that can either influence, or be influenced by, that clade.

All clades can potentially be influenced by temperature and sea level changes. For cetaceans, aside from these two variables, we may expect a dependency on silica levels, archaeplastida, and ostracods. Indeed, the evolution of cetaceans has been influenced by diatoms, which are phytoplankton that require silica [START_REF] Marx | Climate, Critters, and Cetaceans: Cenozoic Drivers of the Evolution of Modern Whales[END_REF]. Archaeplastida include green algae that can have a negative impact on diatom blooming [START_REF] Keating | Blue-Green Algal Inhibition of Diatom Growth: Transition from Mesotrophic to Eutrophic Community Structure[END_REF], and thus potentially also affect cetaceans. Finally, ostracods may also be expected to influence cetacean diversification, either directly as a food source or indirectly as diatom predators. For Portulacaceae, biologically relevant variables are CO 2 , δ 13 C, silica, and Archaeplastida.

Indeed, reductions in CO 2 levels have been linked to radiations of angiosperms during the Cretaceous and the rise of C 4 plants, such as Portulacaceae, during the Miocene [START_REF] Gerhart | Plant responses to low [CO2] of the past: Tansley review[END_REF], which is also reflected in global levels of δ 13 C. Finally, Archaeplastida include land plants that can engage in mutualistic or competitive interactions with Portulacaceae. For ruminants, changes in plant life across the globe have been crucial, such that all variables mentioned for Portulacaceae (CO 2 , δ 13 C, silica, and Archaeplastida) are also relevant through their presumed effect on land plants. Archaeplastida can in addition directly influence diversification in Ruminantia via change in their feeding regime.

In order to illustrate how one could test the relative importance of abiotic versus biotic variables, we selected the best supported abiotic and biotic variables and computed the relative probability of these two models based on their Akaike weights; we also computed the Akaike weight of each model among the set of all models and summed these

Akaike weights over all the abiotic variables (or biotic variables) to obtain an estimate of the overall support for abiotic (or biotic) variables. The first approach has the advantage to not be biased by the total number of variables, or biologically relevant variables, in each category (abiotic and biotic). However, it measures the relative importance of the two most supported abiotic and biotic variables rather than the relative importance of abiotic and biotic variables as a whole. Intermediate approaches could be envisioned where the relative support for abiotic and biotic variables is computed based on a fixed, homogenized number of biologically relevant biotic and abiotic variables.

Additionally, for each phylogeny, we assessed the support of the best-fit paleoenvironmental model with smoothing splines computed with increasingly smaller DOFs. This was done in order to test at what resolution of the paleoenvironmental curve support for environment-dependent diversification started to be lost. We compared all fits using AICc as above.

Results

Performance of temperature-dependent diversification models

Temperature-dependent models were able to accurately recover simulated parameter values. Maximum-likelihood estimates of the parameters of the temperature-dependent λ models were unbiased (Figure 1), even when extinction rates were high (Figure 1D).

Maximum-likelihood estimates of the parameters of the temperature-dependent µ models were likewise unbiased, but they were more variable (Supplemental Figure 1). When both λ and µ were exponentially dependent on temperature, parameter recovery remained reliable (Supplemental Figure 2).

Undersampling did not bias parameter estimates, but affected the uncertainty around the estimates (Supplemental Figure 3). While the median estimates did not deviate significantly from the simulated estimates (P < 0.01), the variance increased with decreasing sampling fraction.

Temperature-dependent speciation models can be correctly distinguished from other models. Each model -constant, time-dependent, and temperature-dependent λ modelswere overwhelmingly recovered when they were the generating model (≥ 75%, Figure 2), even when relative extinction was high (Supplemental Figure 4). Trees simulated with a time-dependency on λ were rarely (∼ 5%) recovered by a temperature-dependent model and vice versa, which is especially significant given the overall tendency for the temperature curve to decrease with time. Heterogeneity in speciation rates across lineages did not lead to a false support for temperature dependency: trees with speciation rate shifts were recovered by constant-rate models 32% of the time, models with increasing speciation rate through time 62% of the time, and temperature-dependent models only 6% of the time. The support for increasing speciation rate models can be explained by the fact that lineages that by chance have high speciation rates leave more descendants, such that the average speciation rate over lineages increases through time. The ability to detect temperature-dependency on µ when it occurred was considerably lower: ∼ 20% of trees simulated with a temperature-dependency on µ were accurately recovered by the generating model (Supplemental Figure 5).

As expected, the ability to recover parameter estimates and the temperature-dependent model improved with tree size. For temperature-dependent λ models, median parameter estimates approximated simulated values precisely as soon as trees had more than 100 tips, and the variance of recovered estimates decreased considerably with tree size (Supplemental Figure 6). With the α λ values considered in these simulations (α = ±0.05), the Akaike weight was well above 0.5 for trees with more than 200 tips, and steadily increased to 0.9 for trees with more than 900 tips (Supplemental Figure 6). For temperature-dependent µ models, median estimates for α µ approximated simulated values for trees with more than 400 tips and the range of estimates decreased considerably for trees with more than 600 tips (Supplemental Figure 7). Median estimates for µ 0 were slightly higher than simulated values for all tree sizes, which is consistent with our previous estimates of µ 0 at α ≥ 0.1 (see Supplemental Figure 1). The Akaike weight of the temperature-dependent µ model increased with tree size, but in agreement with the results above it remained low, below 0.4 even for trees with more than 800 tips (Supplemental Figure 7).

When fitting models with degraded temperature curves, we found that accurate parameter estimates could still be properly recovered, unless the environmental data were very highly degraded (i.e., for a smoothed spline constructed with ≤ 15% of the total DOFs available in generalized cross-validated model (Supplemental Figure 8). Fitting models with a degraded temperature curve slightly reduced the ability to distinguish temperature-dependent models from time-dependent ones, but this ability remained good, with < 15% of the phylogenies simulated under a temperature-dependent model with the most degraded curve (DOFs= 5% of the generalized cross-validated value) finding support for time-dependent rather than temperature-dependent models (Supplemental Figure 8).

Performance of model selection on trees simulated with various environmental dependencies

The ability to detect environment-dependence when it exists does not drastically change with the paleoenvironmental variable considered. For all nine paleoenvironmental variables, environment-dependent speciation can be detected and distinguished from a constant rate and time-dependent effect, as soon as the environmental effect is relatively strong (Figure 3). The ability to distinguish between the simulating environment-dependent model and a pure-birth or time-dependent model increases roughly linearly with α λ , from approximately 20% at α λ = 0.15 to 100% at α λ = 2.2. This result is unlikely to be an effect of tree size, as the median tree size does not vary much across α λ values (185 ± 36). As α λ increases, so does the support for the simulated environmental model, as defined by the Akaike weight support for the environment-dependent, time-dependent, and constant-rate models (Supplemental Figure 9), consistent with expectations. For α λ > 1, curves for biotic palaeodata, which tend to be coarser than abiotic curves, show, on average, higher recovery rates, while for α λ < 1, curves for abiotic palaeodata show, on average, higher recovery rates (Supplemental Figure 10A). There is no such noticeable effect on Akaike weights, though (Supplemental Figure 9). There is no apparent effect of linearity, the level of autocorrelation, the direction or magnitude of the best-fit regression slope, degrees of freedom, or the average rate of change of the environmental curves on our ability to recover how the corresponding environment affected speciation (Supplemental Figure 10B-C).

Environment-dependent extinction is generally not detected, regardless of the strength of the environmental effect: the ability to recover environment-dependent µ is < 10% across all α µ values (Supplemental Figure 11). Even when only trees with more than 400 tips are considered (median tree size = 1183 tips), the ability to recover the correct model does not exceed 20%. Environment-dependent extinction is more often detected as an environmental dependence on λ than on µ (Supplemental Figure 11), but this still happens in less than 20% of the trees.

The ability to distinguish between λ-dependence on different environmental variables is good as soon as the environmental dependence is strong enough: while a considerable percentage of trees were recovered by the wrong environment-dependent model for low α λ (a condition under which the environmental model would likely not be selected when compared to a time-dependent model (Figure 3), the correct model was recovered in the majority of trees for α λ > 1 (Figure 4). For several paleoenvironmental dependencies -specifically, δ 13 C, sea level, and temperature -the simulating model was recovered the majority of the time, even at low α λ values. These three environmental variables were also those that were most often misleadingly detected for trees generated under other environmental dependencies (Figure 4). As above, the median tree size did not vary much across α λ values (249 ± 34).

Testing macroevolutionary hypotheses of environment-dependency on empirical phylogenies

The Cetacean phylogeny was best fit by a pure-birth exponential speciation model with a positive dependency on temperature across all posterior distributions (λ = 0.097e 0.02•T (t) for the consensus phylogeny) (Figure 5). The probability that the abiotic model is the best, as calculated by the Akaike weight for the best-fit abiotic model (temperature) when compared to the best-fit biotic model (ostracod diversity) was 0.72. The sum of Akaike weights across abiotic variables was 0.73 ± 0.12 and across biotic variables 0.27 ± 0.15. We furthermore compared fits of the temperature-dependent λ model with splines smoothed by a range of DOFs (5 -190). Using the maximum DOFs (208), as defined by generalized cross-validation, the temperature-dependent model was significantly better supported than a time-dependent model (∆AICc = 6.9) and other environmental models (∆AICc > 1.3).

The time-dependent model only showed a comparable fit to the tree when the temperature curve was smoothed with DOFs below 10. This suggests that diversification in Cetaceans is best supported by the global trend in the temperature curve, although the decrease in AICc support with reduced DOFs is evidence that diversification is also supported by specific fluctuations within the temperature curve. The Ruminantia phylogeny of ( 9) was best fit by a pure-birth exponential speciation model with a positive dependency on δ 13 C across all posterior distributions (λ = 0.175e 0.7•δ 13 C for the consensus phylogeny) with an

Akaike weight for that model (δ 13 C) when compared to the best-fit biotic model (ostracod diversity) of 0.94 (Figure 6A). The sum of akaike weights across abiotic variables was 0.79 ± 0.06 and across biotic variables 0.11 ± 0.06. The Ruminantia supertree of (12), however, was best supported by a constant-rate birth-death model (λ = 0.11, µ = 1.20e -7) (Supplemental Figure 12). The analysis by ( 12) separated dietary groups and found maximum support for a model with a linear dependency of temperature on λ; in our analysis (which did not separate dietary groups), a model with linear dependency on temperature found considerably poorer support than the constant-rate model (∆AICc = 142.329). Portulacaceae was best fit by a pure-birth exponential λ model with a negative dependency on CO 2 for both the maximum-likelihood and Bayesian phylogeny (λ = 2.98e -0.44•CO 2 for the maximum-likelihood phylogeny) (Figure 6B), with an

Akaike weight for that model (CO 2 ) when compared to the best-fit biotic model (foraminifera diversity) of 0.82. The sum of Akaike weights across abiotic variables was 0.75 ± 0.04 and across biotic variables 0.25 ± 0.04. When the curves of the best-fit models for Ruminantia (using the phylogeny of ( 9)) and Portulacaceae were smoothed with increasingly fewer DOFs and refit to the data, there was little effect on the computed AICc for Ruminantia, but a steep increase in AICc values for Portulacaceae at DOFs below 90%, and then again at DOFs below 80% of total DOFs, suggesting that diversification in Ruminantia is supported by the general trend of the δ 13 C curve, whereas in Portulacaceae it is supported by both the general trend and minor fluctuations of the CO 2 curve. For both phylogenies, AICc support for the best-fit environment-dependent model became poorer than for a time-dependent model only when the DOFs dropped to 2, which effectively made the environment-dependent models equivalent to time-dependent models.

Discussion

The history of Earth has been marked by shifting land masses, volcanic eruptions, and continuous global atmospheric change, as well as a widespread and often dramatic turnover and dispersal of life. There is general recognition that changes to biodiversity over geological time are driven by the environment (26) and species interactions (81). The proportional contribution of various abiotic and biotic drivers within individual clades and across the tree of life, however, is still debated [START_REF] Ezard | The challenges to inferring the regulators of biodiversity in deep time[END_REF]. This is due, in part, to the challenges of estimating the effects of these abiotic and biotic factors on clade diversification. We have tested a phylogenetic approach to address these challenges under a model-based framework that allows us to directly relate speciation and extinction rates in a clade to biotic and abiotic factors in the paleoenvironment. Applying this modeling framework to empirical phylogenies, we identify specific environmental variables that likely played a major role in shaping the evolution of the corresponding clades.

Our simulation results point to three aspects of the use of environment-dependent models that should be considered. Firstly, we show that speciation and extinction parameter estimates do not systematically deviate from true values as undersampling increases (as long as undersampling is accounted for when fitting the models), although the variability of estimates certainly does. Secondly, we find that the accuracy of model selection depends strongly on the strength of the environmental dependence measured by the α λ parameter (effect size) and tree size. Finally, we find notable differences in the ability to recover environmental dependency across distinct paleoenvironmental curves. In particular, in the set of paleoenvironmental curves that we considered, recovering the proper environmental dependence under strong dependence was easier for the biotic variables and temperature in comparison with the other abiotic variables. We could not, however, characterize which features of the time-series data influence its recovery rate. A thorough understanding of how differences in the shape of the paleoenvironmental curves influence our ability to make inferences will require simulating environmental curves with distinct and controlled characteristics.

It has been noted before that, although possible in theory, inferring extinction rates from molecular phylogenies is difficult in practice (see (64; 52; 4; 65)). Inferring the environmental dependence of extinction rates is, not surprisingly, also difficult. Specifically, we find that inferring extinction parameters is very variable for small trees (< 200 tips)

and that accurately distinguishing µ-dependent environmental models from constant-rate or time-dependent models can be unreliable.

The strengths and limits of the environmental models identified here suggest some guidelines for best practice. First, parameter estimates for environment-dependent speciation models start to be accurate for trees with more than 100 tips (although this number, of course, also depends on effect size and the nature of the curve). Second, all results related to environmental dependence in extinction should be treated with care, especially for trees with fewer than 400 tips. Second, environment-dependencies that may not seem biologically relevant can be useful negative controls or used to detect potential hidden variables. In our analyses, we did not find strong support for any biologically irrelevant variable, which gives us confidence in the dependencies for which we did find strong support. Third, if a user is interested in analyzing the effect of an environmental variable not considered here, we recommend that s/he directly analyses the statistical properties of the environment-dependent model applied to that particular variable using simulations. Finally, when testing environmental dependencies in speciation rates, inferred α λ values for the best and second-best fitting models can be consulted; if there is a high false discovery rate for the best-fit model at the α λ value inferred for the second-best-fit model (see Figure 4), then the result should be interpreted cautiously.

A recent work that used the temperature-dependent model [START_REF] Condamine | Macroevolutionary perspectives to environmental change[END_REF] to infer diversification in modern birds found rates of speciation and extinction to be negatively dependent on temperature [START_REF] Claramunt | A new time tree reveals Earth history's imprint on the evolution of modern birds[END_REF]. We can assess how robust this result is in light of our simulation analyses. The temperature curve was smoothed with ∼ 25% of the DOFs inferred by generalized cross-validation. Given our observation that the temperature-dependent model is robust to smoothing, we expect this smoothed curve to be sufficient to generate accurate parameter estimates. The tree has 230 tips, which is sufficient for accurately inferring parameter values for λ, but not necessarily for µ. The result that extinction is temperature-dependent, therefore, should be interpreted with caution. The tree is severely undersampled, with a sampling fraction of ∼ 0.02%, which might produce estimates for both speciation and extinction considerably deviated from their true values. Finally, the temperature-dependent model-fitting was conducted on only one tree (the maximum clade credibility tree), which is not ideal. Dating discrepancies can be sources of potential discordance between results, which is something we account for in the Cetacean and (both) Ruminantia phylogenies by analyzing posteriorly sampled trees.

In agreement with previous macroevolutionary (19; 48) and macroecological [START_REF] Whitehead | Diversity of deep-water cetaceans in relation to temperature: implications for ocean warming[END_REF] work, we find that diversification in Cetaceans was positively associated with temperature.

The best-fit model is a pure-birth model with exponential dependence of λ on temperature.

We show in addition that temperature-dependence is better supported than dependence on other paleoenvironmental variables. It has been suggested that Plio-Pleistocene fluctuations in oceanic temperatures created ecological opportunities for allopatric speciation in delphinids [START_REF] Steeman | Radiation of Extant Cetaceans Driven by Restructuring of the Oceans[END_REF] and this may be why we find diversification in Cetaceans to be dependent on temperature. Notably, the inferred α λ value for the second-best-fit model, ostracod diversity, is 0.27 and, at this α λ value, temperature shows a relatively high false discovery rate (around 40%). Therefore, the temperature-dependent result should be interpreted cautiously. The support we find for a positive dependency of ostracod diversity on diversification in Cetaceans (λ = 0.089 * e 0.193 * ostracoddiversity ) suggests that the role of ostracods as a food source may have positively affected Cetacean speciation.

In Ruminantia, the results depended on the phylogeny we used, which is not surprising given the markedly different datations hypothesized by these two phylogenies.

With the phylogeny of ( 9), we found diversification to be best supported by a pure-birth exponential speciation model with a positive dependency on δ 13 C. This substantiates a major role for C 4 grasses, which became dominant during the Miocene and have been hypothesized to have spurred diversification in ruminant mammals [START_REF] Cantalapiedra | Dietary innovations spurred the diversification of ruminants during the Caenozoic[END_REF]. The δ 13 C data we used were collected from marine samples rather than from terrestrial samples (e.g., bovid tooth enamel), and so are not a direct measurement of the relative abundance of C 4 and C 3 grasses. However, there are notable similarities in these data, particularly during the Neogene [START_REF] Cerling | Global vegetation change through the Miocene/Pliocene boundary[END_REF]. For the Ruminantia supertree of ( 12), we find support for a constant-rate birth-deal model rather than any environment-dependent model, which is contrary to what was found in [START_REF] Cantalapiedra | Dietary innovations spurred the diversification of ruminants during the Caenozoic[END_REF]. This may be because that study conducted separate analyses on different dietary groups.

Our analysis of Portulacaceae, which was best fit by a pure-birth model with a negative exponential dependency on CO 2 , adds statistical evidence to the often described macroevolutionary relationship between CO 2 levels and angiosperm diversification.

Previous work has shown that angiosperm diversity in the Cretaceous was negatively correlated with atmospheric CO 2 [START_REF] Ehleringer | A history of atmospheric CO2 and its effects on plants, animals, and ecosystems[END_REF]. The reason for this negative correlation, however, is debated: some argue that so-called "CO 2 starvation" drove radiations in angiosperms as they competed and searched for new niches (79; 69); others that diversification in angiosperms caused a decline in CO 2 through the geological weathering of calcium and silica (43; 82). Our result cannot resolve this conflict, but does extend the previously observed negative relationship into the Cenozoic for an angiosperm family, and furthermore quantifies the effect of CO 2 levels on that family's diversification. We can be confident that the CO 2 model is not a false positive, as CO 2 is rarely selected when it is not the generating model (Figure 4). In particular, when the second-best-fit model (foraminifera diversity) is the generating model, there is a zero or near-zero selection for the CO 2 model, even at much lower α λ values than the one estimated here (4.47) for the dependency to foraminifera (Figure 4).

For all three phylogenies, we estimated the support of all biotic and abiotic models, regardless of whether they were biologically relevant or not. In the three clades, we found a strong support for the best abiotic model when compared to the best biotic one, as well as stronger support for all abiotic variables combined compared to all biotic variables combined. This result may be biologically true, in which case it emphasizes the overarching importance of the global abiotic environment on species diversification in vastly different organisms. It may instead reflect the lack of relevant biotic variables we tested.

We note that the absence of extinction inferred in Cetacea and Ruminantia is not realistic, and is inconsistent with what we know from the fossil record (63; 12). Accounting for heterogeneity in diversification rates across these groups could help in recovering more meaningful extinction rate estimates [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF]. This could be done by applying the same type of tree partitioning that was used by [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF], while fitting temperature-dependent models.

Even though extinction seems to be poorly estimated in Cetacea and Ruminantia, this does not affect the confidence we can have in our inference of the dependency of speciation rates to environmental changes: our simulations suggest that high levels of extinction do not affect the ability to properly infer environment-dependent speciation, and that environment-dependent extinction is not mistakenly interpreted as environment-dependent speciation.

There are many potential extensions of the environment-dependent models presented here. First, paleoenvironmental variables not included in this work could be considered, such as planetary albedo or the origination rate of marine guilds. Second, models accounting for the simultaneous effect of several paleoenvironmental variables as well as their interaction may provide more resolution on how various abiotic and/or biotic factors combine to influence species diversification. While fitting such models is already possible with the current implementation (any kind of dependence of speciation and extinction rates with environmental variables can be imputed, including dependencies on several environmental variables and their interactions), one would need to test the performance of such models using simulations before trusting the associated empirical results. Implementing the potential influence of hidden variables as well as a variable selection procedure would also be important directions for future work.

The way we have accounted for effects of the biotic environment is by considering interclade biotic interactions (e.g., the effect of ostracod diversity through time). However, there is a wide range of intraclade processes that are thought to influence diversification dynamics. For example, speciation may slow down as ecological opportunity recedes and interspecific competition intensifies following an adaptive radiation, resulting in a diversity-dependent effect (62; 27). Alternatively, speciation may accelerate if intra-clade species interactions (such as competition) promote reproductive isolation through phenotypic divergence (71; 2). Future forms of the environment-dependent model, therefore, could be developed to additively or synergistically account for intraclade and interclade effects. Similar models were recently developed for analyzing fossils (46; 29).

Models directly incorporating intraclade effects on phenotypic divergence have also started to be developed (57; 22). Accommodating such factors in environment-dependent models would be useful for evaluating the relative effect of the environment, intra-and interclade species interactions, as well as their combined effect, on clade diversification.

We can also imagine adapting the model to analyse population genetic models on genealogical trees. Here, the environment-dependent model could be used to analyse how demographics (birth and death rates) have varied with time-dependent environmental variables (e.g., the rate of influx of immigrants in a region, the density of air pollution in a city). In traditional population genetic models, population structure is measured as a function of gene flow (88). Changes in gene flow can be the result of recent increases in population size or rates of dispersal; and so the factors that promote or limit gene flow cannot be easily inferred [START_REF] Slatkin | Gene flow and the geographic structure of natural populations[END_REF]. The environment-dependent model, adapted to analyse the effects of, for example, habitat range contraction, environmental disturbance, or precipitous geographic isolation on birth and death rates within a population may therefore help resolve the contribution of those factors to changes in the population. As above, one would need to test the performance of such models using simulations before trusting the associated empirical results.

phylogenetic comparative toolbox is being extended to better account for the effect of environmental changes. Similar to its use here, where we have shown the effect of different environmental variables on diversification in clades, we have seen a growing body of work relating, for example, decreases in temperature to accelerations in diversification in birds [START_REF] Claramunt | A new time tree reveals Earth history's imprint on the evolution of modern birds[END_REF], increases in sea level to increased extinction in butterflies [START_REF] Condamine | Deciphering the evolution of birdwing butterflies 150 years after Alfred Russel Wallace[END_REF], and even orogeny to Andean plant diversification [START_REF] Pérez-Escobar | Recent origin and rapid speciation of neotropical orchids in the world's richest plant biodiversity hotspot[END_REF]. Models of phenotypic evolution accounting for environmental variations have similarly been developed recently, and have found decreases in temperature to drive increases in the rate of body mass evolution in birds and mammals [START_REF] Clavel | Accelerated body size evolution during cold climatic periods in the Cenozoic[END_REF]. Future analyses with these two types of models should allow us to better understand whether diversification and phenotypic evolution respond similarly (or differently) to past climatic changes, as well as the link between environmental variations, phenotypic divergence, and diversification.

Environment-dependent diversification models provide a statistical forum for testing a range of macroevolutionary hypotheses on species trees, implemented in user-friendly software. They provide a means for identifying which paleoenvironmental variable has most impacted diversification in a clade and how specifically that variable has impacted diversification. We have analyzed the statistical behavior of and provided guidelines for an environment-dependent model of species diversification. We think this will help researchers test fundamental macroevolutionary hypotheses on the effects of various abiotic and biotic factors on species diversification.

λ 0 =0.2 λ 0 =0.25 λ 0 =0.3 μ 0 =0.02 α μ =0.02 μ 0 =0.02 μ 0 =0.06 μ 0 =0.1 λ 0 =0.3 α μ =0.01 α μ =-1 α μ =-0.6 λ 0 =0.5 μ 0 =0.1 α μ =0.4 α μ =0.2 λ 0 μ 0 α μ A B C
Supplemental Figure 1: Recovered parameter estimates for trees simulated with a constant λ and an exponential dependency of µ on temperature. Simulations with: (A) varying λ 0 , constant µ 0 , and constant α µ ; (B) constant λ 0 , varying µ 0 , and constant α µ ; and (C) constant λ 0 , constant µ 0 , and varying α µ . Dashed red lines mark the simulated parameter value.

Akaike weight

Supplemental Figure 4: Akaike weights for constant-rate models, time-dependent models, and temperature-dependent models across all temperature-dependent λ trees simulated with a high extinction rate (λ = 0.75e 0.35t , µ = 0.7). Supplemental Figure 5: The ability to correctly recover trees simulated under constantrate (λ = 0.2, µ = 0.01), positive and negative dependence of µ on time (λ 0 = 0.5, µ 0 = 0.01, α µ = ±0.4), and positive and negative dependence of µ on temperature (λ 0 = 0.4, µ 0 = 0.1, α µ = ±0.4). In addition to the simulating models, trees were fit with models with an exponential dependency of λ on temperature. Each column shows the percentage of trees simulated under the described models that find support.

  3 and α µ = 0.01; and α µ = -1, -0.6, 0.2, 0.4, fixing λ 0 = 0.5 and µ 0 = 0.1. Additionally, we simulated trees with both a λ-and µ-dependency on T, with (1) λ 0 = 0.05, α λ = -0.15, 0.15, µ 0 = 0.05, and α µ = 0.1; (2) λ 0 = 0.25, α λ = 0.4, µ 0 = 0.05, and α µ = 0.1; (3) λ 0 = 0.02, α λ = 0.5, µ 0 = 0.05, and α µ = 0.1; (4) λ 0 = 0.075, α λ = 0.5, µ 0 = 0.05, and α µ = -0.4, -0.15, 0.35, 0.5. Values were chosen to test the most amount of variation in the λ and µ dependency, while reliably simulating surviving trees. For each combination of parameter values, we simulated 5000 trees using the RPANDA function sim_env_bd, stopping the simulations when the equivalent of 52 Myrs was achieved. The median tree size (of trees with more than 50 tips) was 280 for trees with a λ-dependency on T, 435 for trees with a µ-dependency on T, and 812 for trees with both a λ-and µ-dependency on T. Finally, we fitted the generating model (exponential dependence of λ

Figure 1 :Figure 2 :

 12 Figure 1: (A) Recovered parameter estimates for trees simulated with an exponential dependency of λ on temperature. Simulations with (A) varying λ 0 and constant α λ and µ 0 , (B) varying α λ and constant λ 0 and µ 0 , and (C,D) varying µ 0 and constant λ 0 and α λ . Dashed red lines mark the simulated parameter value.

Figure 3 :

 3 Figure 3: Model identifiability for trees simulated with an exponential dependency of λ on different paleoenvironments, X, (λ = λ 0 e αX(t) and µ = 0) for varying values of α λ . (A) Rate of recovery is defined as the percentage of trees with better AICc support for the simulating environment-dependency model versus an exponential time-dependent model and constant rate model. (B) Scaled paleoenvironmental curves. Colors in (A) and (B) correspond. Dashed lines correspond to biotic variables.

Figure 4 :

 4 Figure 4: Rate of model identifiability between environment-dependent models. For each set of trees simulated with an exponential dependency of λ on a particular paleoenvironmental curve, X (λ = λ 0 e αX(t) , µ = 0), the percentage of trees recovered by each paleoenvironmental model for different values of α λ . The simulated model is listed above each barplot. Colors correspond to Figure 3; abiotic variables are circumscribed by the dashed box.

Figure 5 :

 5 Figure 5: Environment-dependency in Cetaceans. (A) AICc support for different environment-dependent models, a constant-rate birth-death model, and an exponential timedependent model (without extinction) on a distribution of 100 posteriorly sampled probabilities of the Cetacean phylogeny. Colors and line type correspond to Figure 3. (A, inset) AICc support for temperature-dependent models fitted with smoothed splines defined by varying DOFs. AICc value for the time-dependent model is shown (dashed line). (B) The consensus Cetacean phylogeny and λ time inferred from the temperature-dependent model (on the consensus tree).

Figure 6 :

 6 Figure 6: Diversification dynamics in (top) Ruminantia and (bottom) Portulacaceae. (A) AICc support for different environment-dependent models, a constant-rate birthdeath model, and an exponential time-dependent model (without extinction). (A, inset) AICc support for the best supported model fitted with smoothed splines defined by varying DOFs. (B) λ over time inferred by the best fit model in (A). (B, inset) λ as a function of the environment for the best fit model, where λ is shown as a function of a non-scaled environmental curve. (C) Consensus phylogeny and (D) representative organisms: (top to bottom) Okapia johnstoni, Cephalophus silvicultor, Gazella gazella, Ovis aries, Kobus leche, Bos taurus, Talinum paniculatum, Portulaca molokiniensis, Pereskia sacharosa, Alluaudia ascendens, Montiopsis andicola. The figures are available here (top to bottom): https://www.flickr.com/photos/dkeats/5269538538/sizes/z/, https://www.flickr.com/photos/7326810@N08/1332824505,https://www.flickr.com/photos/93882360@N07/12801059065,https://commons.wikimed

Supplemental Figure 2 :Supplemental Figure 3 :

 23 Recovered parameter estimates for trees simulated with an exponential dependency of λ and µ on temperature. Dashed red lines mark the simulated parameter value. The effect of undersampling on parameter estimation. Parameter estimates for trees simulated with (A) a positive exponential dependency of λ on temperature and constant extinction (µ = 0.05) and (B) a positive exponential dependency of µ on temperature with constant speciation (λ = 0.5) by fitting the temperature-dependent model. Parameter estimates are shown for trees with increasingly smaller sampling fractions, which were achieved by jackknifing the simulated trees by a fixed % and then fitting the environment-dependent model. Simulated parameters are marked by dashed red lines.

  a t u r e -d e p e n d e n t

Supplemental Figure 6 :

 6 Parameter estimates for trees simulated with a positive (left) and negative (right) exponential dependency of λ on temperature by fitting the temperaturedependent model. Parameter estimates are shown for trees with different species richness. Simulated parameters are marked by dashed red lines. Akaike weights are shown for trees fitted with constant-rate models (light grey), time-dependent models (medium grey), and temperature-dependent models (dark grey). Weights are averaged across all trees within each species richness bracket.

Supplemental Figure 7 :

 7 Parameter estimates for (A) λ 0 , (B) µ 0 , (C) α µ for trees simulated with a positive exponential dependency of µ on temperature by fitting the temperaturedependent model. Parameter estimates are shown for trees with different species richness. Simulated parameters are marked by dashed red lines. (C, inset) A magnified plot of α µ estimates. (D) Akaike weights averaged over all trees for each species richness cohort for models fitted with a constant-rate model (light grey), time-dependent model (medium grey), and model with µ dependent on temperature (dark grey).

Supplemental Figure 8 :

 8 Model selection and parameter estimation in temperaturedependent trees under various spline smoothings in the fitted model. (A) Plots of temperature curve splines smoothed by different degrees of freedom. The time-series of temperature data are shown in grey dots and the smoothed curves in black lines. (B,C) Parameter estimates for trees simulated with an exponential dependency of speciation on temperature (B, λ = 0.2e 0.05•T (t) ; C, λ = 0.3e -0.05•T (t) ) with temperature curves determined using generalized cross-validation (degrees of freedom=208), where the fitted temperature-dependent models have temperature curve splines smoothed by different degrees of freedom. Simulated parameters are marked by dashed red lines. (D) The percentage of temperature-dependent trees, simulated with a temperature curve determined using generalized cross-validation, best supported by models fit with temperature curve splines smoothed by different degrees of freedom versus constant-rate models and time-dependent models with an exponential dependence on speciation.

Supplemental Figure 10 :Supplemental Figure 11 :

 1011 Effect of different characteristics on the rate of recovery of environmental curves. The ability to correctly recover the simulated model at different values of α λ for (A) abiotic and biotic variables (see Figure3) and (B) linear and non-linear environmental curves. (C) Barplots for the slope and intercept for regression models fit to autocorrelation functions for multiple lag-times for each environmental curve; the correlation, slope, and degrees of freedom estimated for generalized least squares (GLS) linear fits to each environmental curve; and the average rate of change of each environmental variable with respect to time. Bars are colored according to panel A. Rate of recovery for trees simulated with an exponential dependency of µ on different paleoenvironments, X, (µ = µ 0 e αµX(t) ) for varying values of α µ and constant λ. (Top) Simulated trees fitted with models with an exponential dependency of µ on the paleoenvironment and a constant λ. (Bottom) Simulated trees fitted with models with an exponential dependency of λ on the paleoenvironment and a constant µ. Colors correspond to Figure3.

Supplemental Figure 12 :

 12 Environment-dependency in Ruminantia computed from the Ruminantia supertree[START_REF] Cantalapiedra | Dietary innovations spurred the diversification of ruminants during the Caenozoic[END_REF]. (A) AICc support for different environment-dependent models, a constant-rate birth-death model, and an exponential time-dependent model (without extinction) on a distribution of 5000 posteriorly sampled probabilities of the Ruminantia supertree. All environment-dependent models have an exponential dependency on the environmental variable. Colors and line type correspond to Figure3.

Acknowledgments We would like to thank Julien Clavel, Jonathan Drury, Guilhem Sommeria Klein, Marc Manceau, Cesar Martinez, and Olivier Missa for helpful comments on the manuscript, Odile Maliet for comments and her code to simulate rate-shifted phylogenies, and Juan Cantalapiedra for supplying the distribution of Ruminantia supertrees. EL would like to thank Evan Charles for helpful discussion. Funding was provided through a European Research Council grant (ERC-CoG-PANDA) attributed to HM.

Akaike weight

Supplemental Figure 9: Mean of the Akaike weights for constant-rate models (light grey), time-dependent models (medium grey), and environment-dependent models (dark grey) across all environment-dependent λ trees with the same α λ (see Figure 3).