Pulsed evolution shaped modern vertebrate body sizes

We developed a maximum-likelihood method for fitting Lévy processes to phylogenetic comparative data using restricted maximum-likelihood estimation (REML), by analyzing the phylogenetically independent contrasts (16) (Materials and Methods). The Lévy processes we apply in this work consist of two components: a Brownian motion and a pure jump process. The Brownian motion is characterized by a rate parameter, ${\sigma }^2$, and the pure jump process is characterized by a Lévy measure, $\nu (\cdot )$, where $\nu (dx)dt$ can be thought of as the probability of a jump with a size in the interval $(x,x+dx)$ occurring in the short time $dt$. Note that this model does not couple pulses of evolution to cladogenesis, as in the classical theory of punctuated equilibrium (17). Instead, pulses may occur at any time, sometimes known as “punctuated anagenesis” (18).

Both Lévy processes with jumps and pure Brownian motion accumulate variance proportionally to time (SI Appendix), leading to speculation that it is impossible to distinguish between pulsed and certain incremental models from comparative data (13, 14). For simulations with moderately sized clades ( $>100$ taxa), we had sufficient power to differentiate pulsed evolution from other Simpsonian modes of evolution. This is due to the impact of rare, large jumps resulting in a heavy-tailed distribution of trait change (SI Appendix). Moreover, we saw low false positive rates for identifying pulsed evolution, even in the presence of phylogenetic error (4% for clades with ∼100 taxa, 7% for clades with ∼300 taxa; SI Appendix).

We assembled comparative datasets from time-scaled tree estimates and body size measurements for 66 clades across five major vertebrate groups (Table 1). We computed the sample size-corrected Akaike information criteria weights (wAICc) (19, 20) for each dataset from a panel of models. We used a Brownian motion (BM) to model incremental phenotypic change, in which a lineage’s phenotype follows a wandering optimum. We used the Ornstein–Uhlenbeck (OU) model to model incremental evolution around a single stationary optimum and the early burst (EB) model to capture adaptive radiation, in which the tempo of evolution slows over time. We also compared two different types of jump processes: the compound Poisson with normally distributed jumps [jump normal (JN)] and the normal inverse Gaussian (NIG). The JN process waits for an exponentially distributed amount of time before jumping to a new value; we use this to represent stasis followed by large-scale shifts between adaptive zones. On the other hand, the NIG process is constantly jumping, with larger jumps requiring longer waiting times; this model captures the dynamics of constant phenotypic change, where the majority of change occurs within the width of an adaptive zone, but rare shifts between adaptive zones are still possible. Finally, we examined the combination of the Brownian and jump processes, BM+JN and BM+NIG. We further modeled intraspecific variation and sampling noise using a normal distribution. Examples of how each different Lévy process models trait evolution can be found in SI Appendix.

Fig. 1 shows that that 32% of clades were best fitted by models of pulsed evolution. BM, EB, and OU were each selected substantially less often (18%, 14%, and 2%, respectively). About 34% of clades lacked decisive support for any single mode of trait evolution. We did not find any evidence that a particular evolutionary mode is enriched in any of the five vertebrate groups ( ${\chi }^2$ test, $P\approx 0.10$). Lévy jump models and Brownian models fitted to any given clade yield nearly identical estimates of process variance (SI Appendix), indicating that the tempo and mode of evolution leave distinct patterns in comparative data.

We speculated that clades would favor models that produced similar patterns in comparative data. To characterize this, we applied a principal components analysis to the vectors of wAICc scores of Fig. 1 and then clustered wAICc profiles using k-means (k = 4). The first three principal components explain 85.6% of the variance. Fig. 2 shows that four clusters form around clades that select models of incremental evolution (BM and OU), explosive evolution (EB), and the two models representing pulsed evolution (JN and NIG). Interestingly, the third component separates clades that favored the two different kinds of jump models we explored, with the most pulsed process (JN) running opposite to the infinitely active processes (NIG). This suggests that with more power it will be possible to choose among different models of pulsed evolution with greater confidence.

Under a microevolutionary model of stabilizing selection, intraspecific variation is distributed normally around the adaptive optimum (21). This intraspecific variation in turn approximates the width of a given adaptive zone (6, 22). When fitting models of trait evolution, we inferred a parameter that corresponds to the SD of intraspecific variation. Because we have only one sample per species, this is a noisy measure of intraspecific variation that also includes measurement error. We therefore regressed our estimates of the intraspecific variation, using bird data with known phenotypic standard deviations, and used the regressed estimates as a measure of intraspecific variation, ${\sigma }_{\text{intra}}$ (Materials and Methods and SI Appendix).

We sought to estimate the rate at which lineages escaped their current adaptive zones via shifts. To do so, we computed the expected waiting time, $\bar{t}(z{\sigma }_{\text{intra}})$, until the lineage jumps at least $z$ intraspecific standard deviations away from the current phenotype (Materials and Methods). Fig. 3A shows these waiting times for the pure Jump Normal process. The value at $z=0$ shows the waiting time for a jump of any size to occur, while we take $z\ge 2$ to be a jump outside of the current adaptive zone, which we define as the interval containing roughly 95% of intraspecific variation per lineage. We found that lineages waited between 1 million and 100 million y before shifting to new adaptive zones ( $z\ge 2$), with a median waiting time of $2.5\times {10}^7$ y.

The waiting time between adaptive shifts depends both on the magnitude of the shift and on the underlying jump model (Fig. 3B). Models with jumps differ in two ways: (i) the relative frequencies of jumps within and between adaptive zones (e.g., $z<2$ vs. $z\ge 2$) and (ii) whether the fitted model is a pure jump process or a jump-diffusion process. JN models experience long periods of stasis, and jumps tend to be large. Thus, small changes within a lineage’s adaptive zone occur on approximately the same timescale as jumps between adaptive zones. In contrast, infinitely active processes, as represented by NIG, jump continuously. Under these models, shifts within an adaptive zone occur frequently, which could reflect rapid adaptation to small-scale environmental perturbations, while larger shifts between adaptive zones occur on much longer timescales. Interestingly, we find that both finite and infinite activity pure jump models predict similar rates of shifting between adaptive zones, indicated by the fact that the solid curves are near each other for $z\ge 2$. When looking at a clade that is best fitted by a jump-diffusion model, such as Anguimorpha (Fig. 3B, Right), we find that they explain most of the morphological disparity among taxa via the BM component of the model, and we see that jumps occur extremely rarely compared with pure jump models. Jump-diffusion models with small diffusion components produce jump size–time curves that are essentially identical to the jump size–time curves for the corresponding pure jump model (compare NIG and BM+NIG for Muroidea).

Because most Lévy processes are known only by their characteristic function, we developed a REML algorithm that operates on characteristic functions. We computed the likelihood of the independent contrasts by proceeding from the tips to the root of the phylogeny. To do so, we recursively update the estimate of the trait at internal nodes as a linear combination of the trait value at the two daughter nodes. We also include an additional term to model the uncertainty in trait values at the tips of the tree due to both intraspecific variation and measurement error. Details of the algorithm can be found in SI Appendix.

To fit the model to data, we used the R programming language and relied on several R packages (37–39). Our package is available at https://github.com/Schraiber/pulsR. Optimal parameter values are obtained by maximizing the likelihood, using the Nelder–Mead simplex method (40). To ensure that we obtained true maximum-likelihood estimates, we performed multiple independent parameter searches (10 for empirical data, 4 for simulated data). In addition, we validated that all maximum-likelihood estimates agreed with model hierarchies; e.g., the maximum likelihood for OU must be greater than or equal to that for BM.

We assessed the power and false positive rate of our method by manipulating tree size and phylogenetic error in a controlled simulation setting. Each dataset was simulated by sampling one tree from an empirical posterior distribution of trees (41) and then simulating trait data under each of the seven candidate models (BM, OU, EB, JN, NIG, BM+JN, BM+NIG). Simulation parameters are given in SI Appendix, Table S1. We then computed AIC weights for each of the seven datasets across the same seven candidate models, once under the originally sampled tree (the “true” tree) and then again under the maximum clade credibility tree that summarizes the posterior (the “consensus” tree). Repeating this procedure 100 times, we counted how often each model type is selected in the presence and absence of tree error. This experiment was conducted for two clades with sizes representative of our empirical datasets (Procellariidae with 106 taxa and Tityridae, Tyrannidae, and allies with 325 taxa). In total, we simulated 1,400 datasets and performed 39,200 maximum-likelihood fittings.

We compiled numerous phylogenetic and trait measurement resources to measure body size evolution (Table 1). The species relationships within each clade were represented using fixed phylogenies with branch lengths measured in millions of years. Body size data were represented by body length for fish, amphibians, and reptiles and, mass being a function of volume, by the cube root of body mass for mammals and birds. Trait measurements were then log transformed, under the assumption that trait evolution operates on a multiplicative scale (42). Handling the data in this manner permits the comparison of evolutionary parameters between clades because our body size and tree data share a common scale. Broadly, we applied the same data assembly strategy to all clades within each of the five vertebrate groups. Clades were generally either extracted from large, taxon-rich phylogenies (amphibians, reptiles, birds) or pooled across independently estimated phylogenies (fish, mammals). Availability and quality of trait data varied between vertebrate groups, each requiring special handling.

We computed sample size-corrected AIC weights (wAICc) for each clade across four model classes: BM, OU, EB, and Lévy processes (Fig. 1). The Lévy process with the highest AICc score was chosen to represent all Lévy processes within the class. Model selection required the model class to be at least twice as probable as any competing model class. With BM being a special case of the remaining six models, this threshold lies on the cusp of the maximum relative probability that it might receive under wAICc. Numerous alternative analyses under various model classifications and assumptions are located in SI Appendix.

For Lévy processes with a jump component, we are interested in the waiting time until a jump outside of the current adaptive zone occurs. Under the symmetric jump models we consider here, the waiting time for a jump larger than $x$ is We characterized the width of adaptive zones by the amount of intraspecific variation in a clade. To provide an accurate estimate of intraspecific variation for all clades, we regressed our observed values of ${\sigma }_{\text{tip}}$ against direct measurements of intraspecific SD in birds (59) ( $n=16$, $P<0.005$, $R^2=0.53$, slope $=0.200$, intercept $=0.059$). Further details for this analysis and results using the uncorrected ${\sigma }_{\text{tip}}$ estimates are given in SI Appendix.