Accelerated body size evolution during cold climatic periods in the Cenozoic

A General Model of Phenotypic Evolution Accounting for Environmental Variations.

To test the effect of past measured environmental variables on rates of phenotypic evolution, we extend the BM process (16⇓⇓–19) with time-varying evolutionary rate σ(t):dX(t)=σ(t)dB(t),[1]where dB(t) is a white noise with mean = 0 and variance of dt. We allow σ(t) to be influenced by one or k environmental variables E1(t),E2(t),…,Ek(t), which themselves vary through time:σ∼(t)=σ(t,E1(t), E2(t), …,Ek(t)).[2]The likelihood corresponding to this model is the classical multivariate normal distribution (18, 40), with the variance–covariance matrix given byVij=∫0Sijσ∼2(t)dt=∫0Sijσ2(t,E1(t), E2(t), …,Ek(t))dt,[3]where Sij represents the time between the root and the most recent common ancestor of species i and j (e.g., ref. 41 has related models). To speed up the computation of the likelihood, we used a stretching–pruning approach, which consists of transforming (stretching) the branches of the tree according to the expected variance–covariance (42) before computing the likelihood recursively using a fast dynamic algorithm based on independent contrasts (pruning) (40, 43, 44). The integrals 3 were computed numerically using the Gauss–Kronrod quadrature formula (45) implemented in the “integrate” function from the stats R base package (46). Finally, maximum likelihood optimization was performed using the quasi-Newton method (47) (L-BFGS-B) implemented in the “optim” function in R. These implementations are available in the RPANDA package (48) publicly available from the CRAN repository (function fit_t_env). They can be used on both ultrametric and nonultrametric trees, therefore allowing the possibility to analyze combined fossil and extant data.

We applied this general model to test if and how rates of phenotypic evolution are related to changes in temperature T. We scaled the temperature curve between zero and one; in what follows, T stands for scaled temperature. We considered two simple models relating phenotypic rates σ2 to temperature T either linearly [σ2(t)=σ02+βT(t) (the Clim-lin model)] or exponentially [σ2(t)=σ02eβT(t) (the Clim-exp model)], where σ02 is the hypothetical clade-specific phenotypic rate at an average global-scale temperature of 0∘C and β reflects the strength and direction of the dependency to temperature. In these models, rates of phenotypic evolution are increasing with increasing temperature when β is positive, are decreasing otherwise, and reduce to a constant rate BM when β= 0. For computational convenience, the Clim-lin model was parameterized as σ2(t)=σ02+(β−σ02)T(t), such that with a scaled temperature curve, σ2 is made up between σ02 and β, increasing with temperature when β>σ02 and decreasing with temperature when β<σ02 (49). We thoroughly tested the ability of our approach to recover input parameters using extensive simulations (SI Appendix).

Model Comparison.

We compared the fit of the climatic models with six competitive models of trait evolution. We fitted the classical BM and an Ornstein–Uhlenbeck (OU) process, which both assume a constant diffusion σ (50, 51). On ultrametric trees and assuming that the root state is at the optimal trait value, the likelihood of the OU model is identical to a time-dependent model with σ increasing exponentially with time [known as the accelerating rate (AC) model (52); we name it the exponential increase (ACexp) model here for clarity]. We, therefore, refer to this process as the OU/ACexp process. We also consider a time-dependent model with σ decreasing exponentially with time [the early burst model (41, 53); we name it the exponential decrease here for clarity] and a time-dependent model with σ varying linearly with time either positively or negatively (4, 41) (coined linear increase or decrease). Finally, we consider the two models that have been used so far to model diversity-dependent effects, with σ constrained to decay with the number of lineages (3, 4) either exponentially (exponential diversity dependence) or linearly (linear diversity dependence). The relative statistical support for the various models was assessed using the Akaike weights (54). We thoroughly assessed the statistical properties of our model comparison framework using intensive simulations (SI Appendix). In particular, we tested our ability to recover the climatic model when it was the generating model and also, that it was not spuriously detected when it was not the generating model.

Body Size Data.

We extracted body mass estimates (in grams) for 9,993 bird species from the EltonTraits 1.0 database (55) and 3,574 mammal species from the PanTHERIA 1.0 database (56). We discarded estimates for 261 bird species that were based on genus or family mean values and that could have biased our evolutionary rate estimates. Body mass estimates were log-transformed before analysis.

Phylogenetic Trees.

Bird phylogenies were taken from the recently updated (v2.iii) (57) posterior distribution by Jetz et al. (58), from which we discarded species that did not have molecular information. Mammal phylogenies were taken from two sources. The first consisted of 1,000 trees sampled from the pseudoposterior distribution by Kuhn et al. (20), which was obtained by randomly resolving polytomies from the widely used supermatrix tree by Bininda-Emonds et al. (59). Because these random polytomy resolutions could inflate evolutionary rate estimates and bias our results (60⇓–62), we also conducted all of our analyses on 1,000 trees from the posterior sample (v.1.002) of a recently published phylogeny of 4,160 extant mammal species by Faurby and Svenning (21) largely based on sequence alignments and ages by Meredith et al. (63). However, Faurby and Svenning (21) focused on resolving topological conflicts rather than branch length, and the authors themselves caution against interpreting branch lengths in their phylogeny. Thus, although the two phylogenies supported consistent results, we reported results for only the birds in the text.

We aligned the phylogenetic and body size data; to test the robustness of our results to taxonomic scale, we conducted analyses at both the order and family levels. We dismissed phylogenies with less than 50 species, because results from simulations showed that a minimum of 50 species was necessary to be able to statistically distinguish our climatic models from other models (SI Appendix). For birds, this alignment resulted in the analysis of 16 orders and 36 families, representing a total of 6,110 species. For mammals, this alignment resulted in the analysis of 12 orders and 15 families, representing a total of 3,355 species [11 orders and 12 families representing a total of 2,664 species for the trees by Faurby and Svenning (21)].

Temperature Data.

We used the temperature curve by Cramer et al. (14). Similar to the more widely used Zachos curve (13), the curve by Cramer et al. (14) is derived from benthic foraminiferal (bf) δ18Obf isotopic ratio. However, contrary to the Zachos curve, the curve by Cramer et al. (14) accounts for fluctuations in sea water (sw) δ18Osw through time, which is important for periods of large-scale glaciations when differences in δ18Osw can go up to −1.11o/oo in Vienna Standard Mean Ocean Water (VSMOW) (64). In addition, the curve by Cramer et al. (14) provides temperature estimates for the last 108 My, thus spanning the full time range over which extant bird and mammal orders originated. Although this curve is derived from the marine record, it correlates well with the more fragmented continental record (65). Rather than local or seasonal fluctuations, these curves reflect planetary-scale climatic trends that are expected to have led to temporally coordinated changes in several clades (7, 9, 66).

Simulating Cope’s Rule.

We simulated Cope’s or Depéret’s rule (23)—the general tendency for increasing body sizes through time—to check whether this trend could artificially favor the support of our climatic model. We simulated evolution toward larger size as taxa chasing an increasing size optimum (67) using a generalization of the OU model (also called Hull–White model) (4):dX(t)=α[θ(t)−X(t)]dt+σdB(t).[4]We simulated an adaptive optima changing either linearly through time according to θ(t)=θ0+μt or linearly as a function of temperature θ(t)=θ0+μT(t) according to Bergman’s rule (15, 22). Our simulations were run on the phylogenies corresponding to each order with two sets of parameter values. We chose α values corresponding to a phylogenetic half-life [time for the OU process to reach one-half the time to stationarity (50)], representing 10 and 100% of the tree height; σ was chosen to be 2α times the observed trait variance [as expected under the stationary condition (50)], and μ was fixed to 0.02 (−0.02 in the case of the optima tracking temperature). Simulations were performed recursively using a forward algorithm from the root to the tips using our own code.

Assessing the Effect of Episodic Bursts of Phenotypic Evolution.

It has been proposed that phenotypic evolution in most vertebrate groups proceeds by rare but substantial bursts along isolated branches (24, 39) [from a 2- to a 52-fold increase, with median value around five in the mammalian supertree (figure 1B in ref. 24)] and that such bursts might drive the support of homogeneous models estimated over entire clades (24). We believe this to be unlikely, because the pattern of interspecific covariances for a process of punctuated evolution with large normally distributed changes is expected to be almost nondifferentiable from that of a BM process (18). We nevertheless tested the possibility that localized shifts in trait evolution artificially favor support for our climatic models using simulations. For each order-level phylogeny, we performed 1,000 simulations of Brownian evolution with localized shifts in σ. Each simulation consisted of randomly selecting edges in the tree where shifts occur (the number of shifts was a proportion ranging from 1 to 10% of the number of species in each order). The amplitude of each shift was drawn from a truncated log-normal distribution with mean = 1.5, variance = 0.5, and lower and upper bounds = 2 and 52, respectively; these parameters reproduce the range of rate increases previously observed on mammals, with a median value around five (24). The simulations were performed by stretching the randomly selected branches according to the selected rate increases before simulating a homogeneous Brownian process with σ2= 1 using the recursive function “rTraitCont” from the R package ape (68). We then fitted eight competitive models to each simulated dataset and compared their relative fits.

Testing Whether Specificities of the Temperature Curve Matter.

To test whether the fit of the climatic model could be explained by the overall cooling trend over the Cenozoic rather than specificities of the temperature curve, we assessed the impact of increasingly smoothing the temperature curve on the support of the Clim-exp model. We used cubic splines with a decreasing effective number of dfs to smooth the curve (69). For each degree of smoothing, we computed the proportion of trees from the posterior distribution of 1,000 trees for which the fit with the smoothed climatic curve remains as good as with the original curve (|AICsmoothed−AICoriginal|< 4) (SI Appendix, Figs. S22–S24). The ΔAIC threshold of four represents a useful approximation for the 95% confidence set on the reference (unsmoothed) model (54).

Assessing the Robustness to ME.

We used simulations to test if ME in the body size data could artificially drive the observed climatic signal. We first derived empirical distributions of ME for each bird order using data from ref. 70 (SI Appendix). Next, for each order, we simulated tip data under BM and OU on 1,000 trees from the posterior distribution, and on each of these tip data, we added ME drawn from the empirical distribution (SI Appendix). Finally, we conducted our model fitting procedure on the resulting simulated data.

To refine our understanding of what type of climatic signal ME would spuriously create or in contrast, blur, we conducted time bin analyses. For each order, we sliced trees from the posterior distribution into 2.5-Ma time bins using the “make.era.map” function in phytools (71). We then jointly estimated maximum likelihood rates for each time bin using the “mvBM” function in mvMORPH (44). We performed these analyses on the empirical body size data and data simulated under BM and OU and with ME as described above. Finally, we reported estimated differences in rates (both empirical and simulated) corresponding to a cold period spanning most of the Oligocene (33.9–25 Ma) and a warm period spanning from the late Oligocene to the mid-Miocene (25–16 Ma). Average rates on these periods were obtained by computing the mean rates across the corresponding time bins (the periods were approximated to span 35–25 and 25–15 Ma, respectively, to match the time bins). This approach is useful to visually inspect temporal trends and focus on specific time periods; however, the uncertainty around estimates in each time bin is high and hampers the statistical assessment of general climatic effects in contrast to our proposed framework.