6.3: Variation in rates of trait evolution across clades - Biology

6.3: Variation in rates of trait evolution across clades - Biology

We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

One assumption of Brownian motion is that the rate of change (σ2) is constant, both through time and across lineages. However, some of the most interesting hypotheses in evolution relate to differences in the rates of character change across clades. For example, key innovations are evolutionary events that open up new areas of niche space to evolving clades (Hunter 1998; reviewed in Alfaro 2013). This new niche space is an ecological opportunity that can then be filled by newly evolved species (Yoder et al. 2010). If this were happening in a clade, we might expect that rates of trait evolution would be elevated following the acquisition of the key innovation (Yoder et al. 2010).

There are several methods that one can use to test for differences in the rate of evolution across clades. First, one can compare the magnitude of independent contrasts across clades; second, one can use model comparison approaches to compare the fit of single- and multiple-rate models to data on trees; and third, one can use a Bayesian approach combined with reversible-jump machinery to try to find the places on the tree where rate shifts have occurred. I will explain each of these methods in turn.

Rate tests using phylogenetic independent contrasts

One of the earliest methods developed to compare rates across clades is to compare the magnitude of independent contrasts calculated in each clade (e.g. Garland 1992). To do this, one first calculates standardized independent contrasts, separating those contrasts that are calculated within each clade of interest. As we noted in Chapter 5, these contrasts have arbitrary sign (positive or negative) but if they are squared, represent independent estimates of the Brownian motion rate parameter (σ2). Basically, when rates of evolution are high, we should see large independent contrasts in that part of the tree (Garland 1992).

In his original description of this approach, Garland (1992) proposed using a statistical test to compare the absolute value of contrasts between clades (or between a single clade and the rest of the phylogenetic tree). In particular, Garland (1992) suggests using a t-test, as long as the absolute value of independent contrasts are approximately normally distributed. However, under a Brownian motion model, the contrasts themselves – but not the absolute values of the contrasts – should be approximately normal, so it is quite likely that absolute values of contrasts will strongly violate the assumptions of a t-test.

In fact, if we try this test on mammal body size, contrasting the two major clades in the tree (carnivores versus non-carnivores, Figure 6.2A), there looks to be a small difference in the absolute value of contrasts (Figure 6.2B). A t-test is not significant (Welch two-sample t-test P = 0.42), but we also can see that the distribution of PIC absolute values is strongly skewed (Figure 6.2C).

There are other simple options that might work better in general. For example, one could also compare the magnitudes of the squared contrasts, although these are also not expected to follow a normal distribution. Alternatively, we can again follow Garland’s (1992) suggestion and use a Mann-Whitney U-test, the nonparametric equivalent of a t-test, on the absolute values of the contrasts. Since Mann-Whitney U tests use ranks instead of values, this approach will not be sensitive to the fact that the absolute values of contrasts are not normal. If the P-value is significant for this test then we have evidence that the rate of evolution is greater in one part of the tree than another.

In the case of mammals, a Mann-Whitney U test also shows no significant differences in rates of evolution between carnivores and other mammals (W = 251, P = 0.70).

Rate tests using maximum likelihood and AIC

One can also carry out rate comparisons using a model-selection framework (O’Meara et al. 2006; Thomas et al. 2006). To do this, we can fit single- and multiple-rate Brownian motion models to a phylogenetic tree, then compare them using a model selection method like AICc. For example, in the example above, we tested whether or not one subclade in the mammal tree (carnivores) has a very different rate of body size evolution than the rest of the clade. We can use an ML-based model selection method to compare the fit of a single-rate model to a model where the evolutionary rate in carnivores is different from the rest of the clade, and use this test evaluate the support for that hypothesis.

This test requires the likelihood for a multi-rate Brownian motion model on a phylogenetic tree. We can derive such an equation following the approach presented in Chapter 4. Recall that the likelihood equations for (constant-rate) Brownian motion use a phylogenetic variance-covariance matrix, C, that is based on the branch lengths and topology of the tree. For single-rate Brownian motion, the elements in C are derived from the branch lengths in the tree. Traits are drawn from a multivariate normal distribution with variance-covariance matrix:

[ extbf{V}_{H_1} = σ^2 extbf{C}_{tree} label{6.5}]

One simple way to fit a multi-rate Brownian motion model is to construct separate C matrices, one for each rate category in the tree. For example, imagine that most of a clade evolves under a Brownian motion model with rate σ12, but one clade in the tree evolves at a different (higher or lower) rate, σ22. One can construct two C matrices: the first matrix, C1, includes branches that evolve under rate σ12, while the second, C2, includes only branches that evolve under rate σ22. Since all branches in the tree are included in one of these two categories, it will be true that Ctree = C1 + C2. For any particular values of these two rates, traits are drawn from a multivariate normal distribution with variance-covariance matrix:

[ extbf{V}_{H_2} = σ_1^2 extbf{C}_1 + σ_2^2 extbf{C}_2 label{6.6} ]

We can now treat this as a model comparison-problem, contrasting H1: traits on the tree evolved under a constant-rate Brownian motion model, with H2: traits on the tree evolved under a multi-rate Brownian motion model. Note that H1 is a special case of H2 when σ12 = σ22; that is, these two models are nested and can be compared using a likelihood ratio test. Of course, one can also compare the two models using AICc.

For the mammal body size example, you might recall our ML single-rate Brownian motion model (σ2 = 0.088, (ar{z}(0) = 4.64), lnL = −78.0, AICc = 160.4). We can compare that to the fit of a model where carnivores get their own rate parameter (σc2) that might differ from that of the rest of the tree (σo2). Fitting that model, we find the following maximum likelihood parameter estimates: (hat{sigma}_c^2 = 0.068), (hat{sigma}_o^2 = 0.01), (hat{ar{z}}(0) = 4.51)). Carnivores do appear to be evolving more rapidly. However, the fit of this model is not substantially better than the single-rate Brownian motion (lnL = −77.6, AICc = 162.3).

There is one complication, which is how to deal with the actual branch along which the rate shift is thought to have occurred. O’Meara et al. (2006) describe “censored” and “noncensored” versions of their test, which differ in whether or not branches where rate shifts actually occur are included in the calculation. In the censored version of the test, O’Meara et al. (2006) omit the branch where we think a shift occurred, while in the noncensored version O’Meara et al. (2006) include that branch in one of the two rate categories (this is what I did in the example above, adding the stem branch of carnivores in the “non-carnivore” category). One could also specify where, exactly, the rate shift occurred along the branch in question, placing part of the branch in each of the two rate categories as appropriate. However, since we typically have little information about what happened on particular branches in a phylogenetic tree, results from these two approaches are not very different – unless, as stated by O’Meara et al. (2006), unusual evolutionary processes have occurred on the branch in question.

A similar approach was described by Thomas et al. (2006) but considers differences across clades to include changes in any of the two parameters of a Brownian motion model (σ2, (ar{z}(0)), or both). Remember that (ar{z}(0)) is the expected mean of species within a clade under a Brownian motion, but also represents the starting value of the trait at time zero. Allowing (ar{z}(0)) to vary across clades effectively allows different clades to have different “starting points” in phenotype space. In the case of comparing a monophyletic subclade to the rest of a tree, Thomas et al.’s (2006) approach is equivalent to the “censored” test described above. However, one drawback to both the Thomas et al. (2006) approach and the “censored” test is that, because clades each have their own mean, we no longer can tie the model that we fit using likelihood to any particular evolutionary process. Mathematically, changing (ar{z}(0)) in a subclade postulates that the trait value changed somehow along the branch leading to that clade, but we do not specify the way that the trait changed – the change could have been gradual or instantaneous, and no amount or pattern of change is more or less likely than anything else. Of course, one can describe evolutionary scenarios that might act like this process - but we begin to lose any potential tie to explicit evolutionary processes.

Rate tests using Bayesian MCMC

It is also possible to carry out this test in a Bayesian MCMC framework. The simplest way to do that would be to fit model H2 above, that traits on the tree evolved under a multi-rate Brownian motion model, in a Bayesian framework. We can then specify prior distributions and sample the three model parameters ((ar{z}(0)), σ12, and σ22) through our MCMC. At the end of our analysis, we will have posterior distributions for the three model parameters. We can test whether rates differ among clades by calculating a posterior distribution for the composite parameter σdiff2 = σ12σ22. The proportion of the posterior distribution for σdiff2 that is positive or negative gives the posterior probability that σ12 is greater or less than σ22, respectively.

Perhaps, though, researchers are unsure of where, exactly, the rate shift might have occurred, and want to incorporate some uncertainty in their analysis. In some cases, rate shifts are thought to be associated with some other discrete character, such as living on land (state 0) or in the water (1). In such cases, one way to proceed is to use stochastic character mapping (see Chapter 8) to map state changes for the discrete character on the tree, and then run an analysis where rates of evolution of the continuous character of interest depend on the mapping of our discrete states. This protocol is described most fully by Revell (2013), who also points out that rate estimates are biased to be more similar when the discrete character evolves quickly.

It is even possible to explore variation in Brownian rates without invoking particular a priori hypotheses about where the rates might change along branches in a tree. These methods rely on reversible-jump MCMC, a Bayesian statistical technique that allows one to consider a large number of models, all with different numbers of parameters, in a single Bayesian analysis. In this case, we consider models where each branch in the tree can potentially have its own Brownian rate parameter. By constraining sets of these rate parameters to be equal to one another, we can specify a huge number of models for rate variation across trees. The reversible-jump machinery, which is beyond the scope of this book, allows us to generate a posterior distribution that spans this large set of models where rates vary along branches in a phylogenetic tree (see Eastman et al. 2011 for details).

Sensory adaptations reshaped intrinsic factors underlying morphological diversification in bats

Morphological evolution may be impacted by both intrinsic (developmental, constructional, physiological) and extrinsic (ecological opportunity and release) factors, but can intrinsic factors be altered by adaptive evolution and, if so, do they constrain or facilitate the subsequent diversification of biological form? Bats underwent deep adaptive divergences in skull shape as they evolved different sensory modes here we investigate the potential impact of this process on two intrinsic factors that underlie morphological variation across organisms, allometry, and modularity.


We use comparative phylogenetic and morphometric approaches to examine patterns of evolutionary allometry and modularity across a 3D geometric morphometric dataset spanning all major bat clades. We show that allometric relationships diverge between echolocators and visually oriented non-echolocators and that the evolution of nasal echolocation reshaped the modularity of the bat cranium.


Shifts in allometry and modularity may have significant consequences on the diversification of anatomical structures, as observed in the bat skull.


Oxygen conductance to the tissues determines aerobic metabolic performance in most eukaryotes but has cost/benefit tradeoffs. Here we examine in lowland populations of a butterfly a genetic polymorphism affecting oxygen conductance via the hypoxia-inducible factor (HIF) pathway, which senses intracellular oxygen and controls the development of oxygen delivery networks. Genetically distinct clades of Glanville fritillary (Melitaea cinxia) across a continental scale maintain, at intermediate frequencies, alleles in a metabolic enzyme (succinate dehydrogenase, SDH) that regulates HIF-1α. One Sdhd allele was associated with reduced SDH activity rate, twofold greater cross-sectional area of tracheoles in flight muscle, and better flight performance. Butterflies with less tracheal development had greater post-flight hypoxia signaling, swollen & disrupted mitochondria, and accelerated aging of flight metabolic performance. Allelic associations with metabolic and aging phenotypes were replicated in samples from different clades. Experimentally elevated succinate in pupae increased the abundance of HIF-1α and expression of genes responsive to HIF activation, including tracheal morphogenesis genes. These results indicate that the hypoxia inducible pathway, even in lowland populations, can be an important axis for genetic variation underlying intraspecific differences in oxygen delivery, physiological performance, and life history.

Figure S1. Predicted mRNA secondary structure and accessibility of the miR-71 target site in the region of the Sdhd transcript flanking the M and I indel alleles.

Figure S2. Minority allele frequencies in contigs of Finland-derived sequences in M. cinxia transcriptome assembly 2.0 (Wheat et al. 2011).

Table S1. Summary of adult butterflies used for experiments associating Sdhd alleles with phenotypes.

Table S2. Sdhd genotype frequencies in samples from three distinct clades from Finland, France, and Spain.

Table S3. Analyses of Sdhd transcript abundance (N = 113) and SDH enzyme rate (N = 50).

Table S4. SDH enzyme activity in relation to 3′-UTR indel genotype a nonsynonymous polymorphism at amino acid 72.

Table S5. Factors affecting metabolic rate during consecutive 10-min episodes of continuously stimulated flight in two different gas mixtures (21% and 14% O2).

Table S6. Mitochondrial integrity in relation to physiological variables and Sdhd alleles.

Table S7. Factors affecting metabolic rate during 3-min episodes of continuously stimulated flight in two different gas mixtures (21% and 14% O2) presented in random order on consecutive days in a sample of 30 male and female butterflies from Spain.

Table S8. Factors related to flight endurance (total CO2 emitted) during the second flight in respirometry experiment 2.

Table S9. Tracheal morphogensis genes differentially expressed (P < 0.05) in Sdhd M vs. non-M butterflies (data from Wheat et al. 2011).

Table S10. Sequences for amplicons amplified from the cDNA of Melitaea cinxia butterflies in quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) assays.

This material is available as part of the online article from: (This link will take you to the article abstract).

Filename Description
EVO_12004_sm_suppmat.zip792.8 KB Supporting info item

Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.

Ecophylogenetics Clarifies the Evolutionary Association between Mammals and Their Gut Microbiota

Our knowledge of how the gut microbiome relates to mammalian evolution benefits from the identification of gut microbial taxa that are unexpectedly prevalent or unexpectedly conserved across mammals. Such taxa enable experimental determination of the traits needed for such microbes to succeed as gut generalists, as well as those traits that impact mammalian fitness. However, the punctuated resolution of microbial taxonomy may limit our ability to detect conserved gut microbes, especially in cases in which broadly related microbial lineages possess shared traits that drive their apparent ubiquity across mammals. To advance the discovery of conserved mammalian gut microbes, we developed a novel ecophylogenetic approach to taxonomy that groups microbes into taxonomic units based on their shared ancestry and their common distribution across mammals. Applying this approach to previously generated gut microbiome data uncovered monophyletic clades of gut bacteria that are conserved across mammals. It also resolved microbial clades exclusive to and conserved among particular mammalian lineages. Conserved clades often manifest phylogenetic patterns, such as cophylogeny with their host, that indicate that they are subject to selective processes, such as host filtering. Moreover, this analysis identified variation in the rate at which mammals acquire or lose conserved microbial clades and resolved a human-accelerated loss of conserved clades. Collectively, the data from this study reveal mammalian gut microbiota that possess traits linked to mammalian phylogeny, point to the existence of a core set of microbes that comprise the mammalian gut microbiome, and clarify potential evolutionary or ecologic mechanisms driving the gut microbiome's diversification throughout mammalian evolution.IMPORTANCE Our understanding of mammalian evolution has become microbiome-aware. While emerging research links mammalian biodiversity and the gut microbiome, we lack insight into which microbes potentially impact mammalian evolution. Microbes common to diverse mammalian species may be strong candidates, as their absence in the gut may affect how the microbiome functionally contributes to mammalian physiology to adversely affect fitness. Identifying such conserved gut microbes is thus important to ultimately assessing the microbiome's potential role in mammalian evolution. To advance their discovery, we developed an approach that identifies ancestrally related groups of microbes that distribute across mammals in a way that indicates their collective conservation. These conserved clades are presumed to have evolved a trait in their ancestor that matters to their distribution across mammals and which has been retained among clade members. We found not only that such clades do exist among mammals but also that they appear to be subject to natural selection and characterize human evolution.

Keywords: Gut microbiome bioinformatics ecology evolution phylogeny taxonomy.

Copyright © 2018 Gaulke et al.


An ecophylogenetic approach to taxonomy…

An ecophylogenetic approach to taxonomy can discover ecologically relevant units of microbial taxa.…

The phylogenetic distribution of conserved…

The phylogenetic distribution of conserved bacterial clades reveals associations between gut microbiota and…

A codiversifying clade within the…

A codiversifying clade within the Bacteroidales contains subtending clades that are unique to…

Primate gut microbiomes have diversified…

Primate gut microbiomes have diversified in a manner correlated with their evolutionary history.…


This research was supported by grants from the Department of Energy’s Biological Systems Science Division (No. DE-SC0016207), Program in Genomic Science and the National Science Foundation (Nos. DEB-1645596 and DEB-1241094). Work at Lawrence Livermore National Laboratory (LLNL) was funded by the Department of Energy through the Genome Sciences Program under contract Nos. SCW1424 and SCW1590, and performed under the auspices of LLNL under Contract No. DE-AC52-07NA27344.


Insects contain ∼ 62% of all ∼ 1.6 million described living species 1 , but the causes of the exceptional diversity of this relatively recent branch ( ∼ 500 million years old) on the Tree of Life remain highly uncertain. Many authors have suggested that herbivory is a major driver of insect diversification 2,3,4,5,6,7,8 , especially feeding on living tissues of vascular plants (although some have also questioned its importance 9,10 ). Thus, a strong relationship between herbivory and insect diversification, if supported, may be key to understanding the overall biodiversity of life on the Earth. The increasing availability of insect phylogenies (both within and between orders 10,11,12,13 ) provides exciting new opportunities to test this relationship.

Given the potential importance of the herbivory-diversification relationship in insects, quantitative tests have been surprisingly few. In their classic study, Mitter et al. 5 compared the species richness of 13 predominantly herbivorous clades with their non-herbivorous sister groups. Most comparisons were between sister clades (for example, families) within orders, including beetles. They found a significant positive impact of herbivory on diversification, based on the higher richness of the predominantly herbivorous clade in most sister-clade comparisons. Other studies have analysed the impact of transitions to phytophagy on diversification rates within beetles (Coleoptera the most species-rich order 14 ). These studies used a similar approach (comparing richness of sister clades). The first supported the role of herbivory in driving beetle diversification 6 , but a later, larger-scale study did not 10 . A paleontological analysis also questioned whether herbivory drives insect diversification, but based on the number of insect families over time rather than on patterns of species richness or diversification 9 . More recent studies have emphasized wings and complete metamorphosis (holometaboly) as drivers of insect diversification 12,15 , but without quantitative comparison of their effects relative to those of herbivory. Thus, the importance of herbivory to broad-scale insect diversification is currently uncertain, despite extensive work relating herbivory (and shifts in host plant species) to speciation at smaller phylogenetic scales 16,17,18 .

Here, we take advantage of increasing phylogenetic information for insects 10,11,12,13 and more powerful comparative methods to address the relationship between herbivory and insect diversification across scales. Our results support the importance of herbivory in increasing insect diversification and helping to shape patterns of species richness among insect orders. However, the impact of herbivory on diversification is more variable within orders, with significant effects within some clades (Diptera, Hemiptera) but not others (Coleoptera, Hymenoptera, Orthoptera). Sister-clade comparisons support these two basic patterns. Wings and holometaboly may also be important drivers of insect diversification, but unlike herbivory, their impact depends strongly on the phylogeny considered. Thus, we demonstrate the importance of herbivory for insect diversification, while also showing that it is not the sole explanation for the remarkable diversity of insects and arthropods. More broadly, our results show that local-scale ecology (e.g., diet, species interactions) helps drive large-scale patterns of species richness in the Earth’s most diverse clade of organisms and over a time span of ∼ 500 million years.

Phylogenetic patterns of trait and trait plasticity evolution: Insights from amphibian embryos

Environmental variation favors the evolution of phenotypic plasticity. For many species, we understand the costs and benefits of different phenotypes, but we lack a broad understanding of how plastic traits evolve across large clades. Using identical experiments conducted across North America, we examined prey responses to predator cues. We quantified five life-history traits and the magnitude of their plasticity for 23 amphibian species/populations (spanning three families and five genera) when exposed to no cues, crushed-egg cues, and predatory crayfish cues. Embryonic responses varied considerably among species and phylogenetic signal was common among the traits, whereas phylogenetic signal was rare for trait plasticities. Among trait-evolution models, the Ornstein-Uhlenbeck (OU) model provided the best fit or was essentially tied with Brownian motion. Using the best fitting model, evolutionary rates for plasticities were higher than traits for three life-history traits and lower for two. These data suggest that the evolution of life-history traits in amphibian embryos is more constrained by a species' position in the phylogeny than is the evolution of life history plasticities. The fact that an OU model of trait evolution was often a good fit to patterns of trait variation may indicate adaptive optima for traits and their plasticities.

Keywords: Anaxyrus Hyla Lithobates Pseudacris Rana phylogenetic inertia.

© 2018 The Author(s). Evolution © 2018 The Society for the Study of Evolution.


Biodiversity Data Assembly.

Using a phylogeny for Saxifragales inferred from a dataset of 301 genes and 627 species covering all major lineages (SI Appendix, Fig. S1) as a backbone constraint for a supermatrix analysis, we generated a tree comprising 1,736 species (72% of species diversity for Saxifragales SI Appendix, Fig. S2), of which 1,455 (61% SI Appendix, Fig. S3) were matched with occurrence data (736,703 occurrence records, an average of 317 occurrences per species). Major relationships recovered were consistent with recent phylogenetic work (31), including monophyly for all clades recognized as families. For the 23 phenotypic traits collected, 1,388 species (58%) had the minimum of 20% trait coverage we imposed for downstream analysis.

Niche and Phenotype Ordination.

The primary loadings on ordinated niche data captured temperature (SI Appendix, Table S1). Low values (Fig. 1, red and yellow) are associated with hotter and, to some extent, wetter habitats, as well as the hottest arid habitats high values (Fig. 1, blue and green) are associated with cooler and drier habitats. Ordinated phenotype data (SI Appendix, Table S2) primarily captured a complex set of traits including plant height woodiness inflorescence, sepal, and petal shape descriptors seed length and several meristic traits (stamen, petal, and sepal number). Low ordinated phenotype values captured primarily low herbaceous rosette plants, and high values captured trees and shrubs.

Ancestral reconstruction across Saxifragales for PC1 of our dataset of 35 environmental variables branches are colored in a rainbow scale from low ordinated values (red and yellow hotter and to some extent wetter habitats, as well as the hottest arid habitats) to high ordinated values (green and blue mostly colder and drier habitats). Black dots at nodes represent major ecological niche shifts (the upper 95th percentile of node–parent node differences). Red dots at nodes represent diversification shifts in the maximum credibility set. The inset density curves show tip rates for diversification (green), niche (orange), and phenotype (blue), all scaled from the minimum to the maximum reconstructed value. Around the edge are photographs of major representative habitats. Family codes are as follows: (a) Peridiscaceae, (b) Paeoniaceae, (c) Daphniphyllaceae, (d) Cercidiphyllaceae, (e) Altingiaceae, (f) Hamamelidaceae, (g) Iteaceae, (h) Grossulariaceae, (i) Saxifragaceae, (j) Cynomoriaceae, (k) Tetracarpaeaceae, (l) Aphanopetalaceae, (m) Penthoraceae, (n) Haloragaceae, and (o) Crassulaceae.

Niche and Phenotype: Conservatism and Correlated Evolution.

Although niche shifts are numerous throughout Saxifragales (Fig. 1), we found strong evidence for phylogenetic constraints on niche space, consistent with some degree of conservatism of both niche and phenotypic traits over >100 My of evolutionary time (λ test niche: P = 8.3e-222 phenotype: p ∼ 0 SI Appendix, Figs. S4 and S5). To discern whether the phenotypic data captured potential niche-delimiting phenotypes, we asked whether the niche and phenotypic predictors were tightly correlated we found strong evidence for correlated evolution between niche and phenotype (P = 4.1e-36) and their pairwise disparity (P < 2.2e-16 SI Appendix, Fig. S6).

Timing of Macroevolutionary Rates.

Net species diversification and rates of niche and phenotypic evolution (lability) as reconstructed by BAMM (Bayesian Analysis of Macroevolutionary Mixtures) all experienced strong increases toward the present, beginning in approximately the mid-Miocene (Fig. 2A and SI Appendix, Fig. S7), contemporaneous with the mid-Miocene Climatic Optimum (32 ⇓ –34). However, surprisingly, diversification rates through time show significantly earlier increases than for phenotypic and niche lability. Comparing the time to 50% of contemporaneous rates, we found diversification, niche, and phenotypic rates had significantly different timing [Fig. 2B ANOVA P < 2e-16 all pairwise rate differences significant for Tukey honestly significant difference (HSD)], with niche and phenotypic rates lagging behind diversification. Phenotype also experienced a smaller lag either before, or more often after, niche rates for this and other analyses, although this ordering was sensitive to two alternative time calibration methods (Methods) and more often niche and phenotype timings were indistinguishable. Consistent with our rate estimates through time, shifts in diversification regimes were generally phylogenetically deeper than major ancestral shifts (upper 95th percentile) in niche and phenotype space (Fig. 2C ANOVA P = 0.00036 all pairwise differences significant for Tukey HSD except niche vs. phenotype). Most major niche shifts (88.1% of the upper 95th percentile) and major phenotype shifts (84.6% of the upper 95th percentile) postdated 15 Mya.

(A) Median rates and rate distributions for net diversification, niche, and phenotypic macroevolutionary rates (colors shown in legend). The unit of diversification is speciation events per million years niche and phenotypic rates are unitless. For relative comparability, the y axis is scaled from zero to the maximum median rate for all datasets. The gray curve in the background is a global temperature dataset (37). High niche evolutionary rates result from large scaling of the PCA ordination of environmental data (SI Appendix, Fig. S41). (B) Box plot showing the distributions of times to 50% of contemporaneous evolutionary rates from the mid-Miocene (15 Mya to present). (C) Box plot showing the distributions of times to present for either major shifts in ancestral reconstructions (environment, phenotype 95th percentile of node–parent node differences) or shifts in the best shift configuration (diversification). (D) Box plot of evolutionary rates for MS (diversification) and σ 2 (niche and phenotype) as proportionally scaled to the maximum (contemporary) rates. Higher clade values for, for example, diversification indicate earlier rate increases. MS here parameterizes extinction as ε = 0.5 see SI Appendix, Figs. S8–S11 for evaluation of extinction fractions.

We also calculated the MS diversification statistic (35) and Ornstein–Uhlenbeck (OU) σ 2 trait rate parameters (36) for the entire tree and every subclade. In addition to the pure-birth MS statistic (ε = 0), we used a series of extinction proportions (ε = 0.1, 0.5, 0.9) to assess the impact of extinction on the MS results. We identified lags in evolutionary rates with early diversification increases followed by niche and then phenotype (Fig. 2D and SI Appendix, Figs. S8–S13), which was significant for all pairwise comparisons (Tukey HSD rates scaled to a proportion of their contemporary rate for comparability P < 0.002 for all comparisons). Higher extinction proportions reduced overall diversification rate scaling but pairwise comparisons remained significant.

Correlation Between Macroevolutionary Rates and Climate/Climatic Variation.

After the mid-Miocene Climatic Optimum, the timeframe in which we observed major increases in macroevolutionary rates (for diversification, niche, and phenotype), there was a very strong correlation between historical global temperature data (37) and macroevolutionary rates reconstructed by BAMM (fit using an exponential model net diversification: F test P < 2.2e-16 R 2 = 0.8352 niche lability: P < 2.2e-16, R 2 = 0.845 phenotypic lability: P < 2.2e-16, R 2 = 0.779). We also observed a weaker, but significant, relationship between rates and variation in temperature (measured as the SD over moving 0.1-My intervals, fit using an exponential model net diversification: F test P < 2.2e-16, R 2 = 0.2003 niche lability: P = 3.23e-2, R 2 = 0.03233 phenotypic lability: P = 6.25e-8, R 2 = 0.05527 SI Appendix, Figs. S14–S19). A combined model incorporating both temperature and temperature variability best explained the data (mean adjusted R 2 increase 0.0163, Akaike weight 0.705 for diversification, ∼1 for niche and phenotype), suggesting that temperature and temperature variability have independent explanatory power.

We fit likelihood models of temperature (n = 9) and time dependence (n = 9) in RPANDA [refs. 38 ⇓ –40 see Fig. 4 and SI Appendix, Tables S3–S5] to evaluate the association of diversification, niche, and phenotypic rates with historical temperature. For niche and phenotypic rates, support for temperature dependence was decisive (combined Akaike weights of linear and exponential models ∼1) Akaike weights moderately supported an exponential temperature relationship with niche lability (0.84), but support for either linear or exponential dependence was equivocal for phenotypic lability. For diversification rates, RPANDA was unable to distinguish a single best model. The best model (exponential speciation constant extinction with respect to time) did not include climate, but ΔAIC (Akaike information criterion) was only 0.325 and the ΔAkaike weight 0.025 compared with the second-best model (linear speciation and extinction with respect to temperature). The sum of Akaike weights for models with temperature was 0.6837 vs. 0.3163 for those without it (Fig. 4 and SI Appendix, Tables S3–S5 and S8–S10), indicating that the majority of the likelihood was in temperature-dependent models.

Robustness to Priors.

Concern has been raised about the sensitivity of BAMM to prior specification (41) to assess this issue, we explored a series of extreme prior formulations on expected rate events and rate priors, varied over two orders of magnitude compared with the recommended BAMM priors (for further details see SI Appendix, Supporting Information Text and Figs. S20–S22). Expected event priors did not produce noticeable effects on rate curves. Rate priors had an effect particularly on early diversification estimates >40 My before present, with higher rate priors especially for niche and phenotypic analyses tending to result in an early burst of evolution ∼80 to 100 My before present that was absent with smaller prior values. This area of the rate curves was sensitive to dating procedure as well as prior (Methods and SI Appendix, Fig. S23), in sum suggesting greater uncertainty deeper in time. However, within the period of interest for this work (≤15 My before present), rate curves were essentially identical to those observed in recommended priors, and overall shape was robust to prior formulation and dating approach.

Niche- and Phenotype-Associated Diversification.

To assess whether observed rate patterns were associated with spatial variation in climate, we ran a series of trait-associated diversification tests with summary geographic data, using both the BAMM-based STRAPP (STructured Rate Permutations on Phylogenies) statistic and the semiparametric es-SIM statistic. We did not see significant associations of elevation (STRAPP: P = 0.868 es-SIM: P = 0.849) or latitude (STRAPP: P = 0.818 es-SIM: P = 0.194) with net diversification. Using STRAPP, we also did not find relationships of elevation and latitude with niche lability (elevation: P = 0.816 latitude: P = 0.101) or phenotypic lability (elevation: P = 0.386 latitude: P = 0.193). We also did not see a relationship between continental biogeography and diversification rates (STRAPP: P = 0.879), niche lability (STRAPP: P = 0.056), or phenotypic lability (STRAPP: P = 0.751).

These results are consistent with a worldwide response of Saxifragales to climatic cooling and aridification, without a clear association with temperate biomes alone (but see below). Ancestral reconstructions of climate space suggest that Saxifragales were already present in temperate habitats before Miocene cooling the majority of the shifts into the most extreme areas of niche space (the upper 97.5th and lower 2.5th percentiles of contemporary values representing, respectively, the hottest desert and equatorial habitats and the coldest polar and montane habitats) are within the last 5 My (niche: 56.72%, phenotype: 61.54% Fig. 1 and SI Appendix, Figs. S24 and S25).

Subclade Patterns.

We generated a series of subclade rate plots for all families with more than 15 sampled species (Fig. 3A) to assess whether the global pattern we observed was general or driven by particular subclades. Diversification timing was equivocal for Hamamelidaceae and Grossulariaceae, with similar curves to niche rates after 15 My, but other examples clearly show that diversification rates predate shifts in other rates, with the timing consistent with that reported for the whole tree. Likewise, while niche and phenotype have a diverse set of timing patterns, patterns in each of two large clades (Saxifragaceae alliance and Crassulaceae alliance cf. Fig. 3B where these taxa are defined, together comprising 97% of species diversity) show a precedence of niche rates consistent with our global analyses. Subclade plots also suggest a negative relationship between niche lability and clade median temperature (measured as the most recent common ancestor ancestral reconstruction of mean annual temperature SI Appendix, Fig. S26). Consistent with this finding, we found that niche lability had a moderately significant association with habitat (STRAPP: P = 0.033). Other cladewise plots did not reveal a clear pattern for net diversification and phenotypic lability (SI Appendix, Figs. S27 and S28).

(A) Clade rates for major clades (all families with greater than 15 sampled taxa). Rate colors are shown in the legend the yellow bars give the date of 15 My representing the Mid-Miocene Climatic Optimum. (B) Box plots showing the distribution of rate shifts for two large subclades of similar size, Crassulaceae alliance (Aphanopetalaceae + Crassulaceae + Haloragaceae + Penthoraceae + Tetracarpaeaceae) and the Saxifragaceae alliance (Grossulariaceae + Iteaceae + Saxifragaceae). For phylogenetic distribution of these shifts see SI Appendix, Supporting Information Text and Figs. S29–S31 and S37–S39.


2.1 Tip rate metrics

where is the tip rate for species i, Ni is the number of branches between species i and the root, bj is the length of branch j, starting at the terminal branch (j = 1) and ending with the root. Jetz et al. ( 2012 ) demonstrated that, for trees deriving from a Yule process, and with mild extinction, the mean λDR across tips converges on the true speciation rate.

We also considered a simpler metric, node density (Freckleton et al., 2008 denoted by λND). This is simply the number of splitting events along the path between the root and tip of a phylogeny, divided by the age of the phylogeny. While λDR down-weights the contribution of branch lengths that are closer to the root, λND equally weights the contributions of all branches along a particular root-to-tip path, regardless of where they occur in time. Under a pure-birth model (μ = 0), both λDR and λND should yield unbiased estimates of the rate of speciation.

The third measure we considered is the inverse of the terminal branch lengths (λTB). Rapid speciation rates near the present should be associated with proportionately shorter terminal branches smaller values of λTB should thus characterize species with faster rates of speciation. This measure has recently been used as a summary statistic to assess model adequacy in trait-dependent diversification studies (Bromham et al., 2016 Gomes et al., 2016 Harvey & Rabosky, 2017 ). Following Steel and Mooers ( 2010 ), we note that the terminal branch lengths can be used to derive an estimate of the speciation rate this follows from the fact that interior and terminal branches have the same expected value under the Yule process (Steel & Mooers, 2010 ). The corresponding estimator for the i’th tip, λTB is approximately 1/2b where b is the length of a given terminal branch (Steel & Mooers, 2010 ). To our knowledge, λTB has not been used to explicitly estimate tip rates as we do here, but given its utility as a summary statistic and general theoretical properties (Steel & Mooers, 2010 ), we see value in comparing the performance of this metric to others currently in use.

Finally, we considered a Bayesian, model-based approach to estimating tip rates. BAMM (Rabosky, 2014 ) assumes that phylogenies are generated by a set of discrete diversification regimes. Using MCMC, the program simulates a posterior distribution of rate shift regimes, from which marginal posterior rate distributions can be extracted for each tip in the phylogeny. Priors for BAMM analyses were set using default settings from the setBAMMpriors function from BAM Mtools (Rabosky, Grundler, et al., 2014 ). The prior parameterizations specified by this function ensure that the prior density on relative rate changes across the tree is invariant to the scale of the tree (e.g. multiplying branch lengths by 10 6 will not change inferences about relative rates across the tree). We denote BAMM tip speciation rates (mean of the marginal posterior) as λBAMM. As BAMM also estimates extinction rates for each regime, we also calculated tip-specific net diversification rate as λBAMMμBAMM, denoted as rBAMM.

2.2 Tip rate metrics estimate speciation, not net diversification

As suggested previously (Belmaker & Jetz, 2015 supplemental analyses in Jetz et al., 2012 ), DR and presumably other tip-based measurements, more accurately estimate the rate of speciation than the rate of net diversification. However, numerous studies continue to refer to DR as a measure of net diversification (Marin & Hedges, 2016 Oliveira et al., 2016 Cai et al., 2017 Quintero & Jetz, 2018 and many others). This is incorrect and it is straightforward to demonstrate that λTB, λND and λDR are more reliable measures of speciation rates and not net diversification rates, at least when extinction is moderate to high.

To illustrate this property of the metrics, we applied all approaches to constant-rate birth–death phylogenies simulated across a range of extinction fractions (ε = μ/λ), including pure-birth trees (ε = 0) as well as trees exhibiting very high turnover (ε = 1). To evaluate accuracy of speciation estimates as a function of ε, we generated 1,000 phylogenies with 100 tips each, where λ and ε were drawn from uniform distributions (λ: [0.05, 0.3] ε: [0, 1]). Importantly, when λ is sampled uniformly with respect to ε, the distribution of r is not uniform: the mean, range and variance in r decrease dramatically as ε increases. To evaluate the accuracy of r as a function of ε, we thus generated a second set of trees by sampling r and ε from uniform distributions (r: [0.05, 0.3], ε [0, 1]). As a result, λ has constant mean and variance with respect to ε in the first set of simulations, and the same is true for r in the second set of simulations (Figure S1). All phylogeny simulations were conducted with the TreeSim package in r (Stadler, 2011 ).

2.3 Assessment of tip rate metrics

We tested the performance of the metrics by compiling publicly available datasets from a number of simulation-based studies (Table 1). By focusing on simulations from previously published work, we thus ensured that the simulation process itself was effectively blinded to the objectives of this study. We further note that our trial datasets included several studies that were critical of BAMM (Meyer & Wiens, 2017 Moore, Höhna, May, Rannala, & Huelsenbeck, 2016 ). These simulated trees include rate heterogeneity in time and across lineages. Together, these phylogenies present a wide range of tree sizes and diversification rate shifts, providing an ideal comparative dataset for our purposes. To more easily distinguish between these tree types in the text, we refer to the BAMM-type, multi-regime time-constant phylogenies simply as ‘multi-regime’, and the multi-regime diversity-dependent phylogenies simply as ‘diversity-dependent’, even though discrete rate shifts are present in both types of trees. In addition to discrete-shift scenarios (e.g. BAMM-type process), we simulated phylogenies under an ‘evolving rates’ model of diversification (Rabosky, 2010 as corrected in Beaulieu & O'Meara, 2015 ) to explore performance of tip rate metrics when diversification rates change continuously and independently along branches, as might occur if diversification rates are correlated with an underlying continuous trait (FitzJohn, 2010 ). In these simulations, we allowed the logarithm of λ to evolve across the tree under a Brownian motion process, while holding ε constant. The magnitude of rate heterogeneity among branches is controlled by the diffusion parameter σ, where greater values lead to greater heterogeneity in speciation rates. Although published phylogenies with rate data were unavailable for this simulation scenario, we used simulation code and parameters taken directly from Beaulieu and O'Meara ( 2015 ) to generate trees with similar statistical properties to those in their study. Simulations were performed with the following parameters: λ = 0.078, 0.103, 0.145, 0.249 and ε = 0.0, 0.25, 0.50, 0.75. We simulated 100 phylogenies for each (λ, ε) pair, and for three values of σ (σ = 0.03, 0.06, 0.12). We evaluated tip rate accuracy by comparing estimated to true tip rates, using the absolute and proportional error metrics described above. We also examined the correlation between true and estimated tip rates, combining tip rates from all phylogenies generated under the same class of diversification process, and visualizing these data as density scatterplots, generated with the LSD package in R (Schwalb et al., 2018 , where colours indicate the density of points.

Simulation model Number of trees Tree size Regime number Source
Single-regime, constant-rate birth–death 100 100 1 Mitchell and Rabosky ( 2016 )
Single- and multi-regime, constant-rate birth–death 100 51–148 1–6 Moore et al. ( 2016 )
Single- and multi-regime, constant-rate birth–death 400 10–4,296 1–67 Rabosky et al. ( 2017 )
Multi-regime, constant-rate birth–death 20 939–3,708 11 Meyer & Wiens ( 2017 )
Single- and multi-regime, constant-rate birth–death 188 4–3,955 1–73 Mitchell et al. ( 2018 )
Single-regime, constant-rate birth–death, lambda uniform 1,000 100 1 This study
Single-regime, constant-rate birth–death, net diversification uniform 1,000 100 1 This study
Pure birth root regime, 1–4 discrete shifts to diversity-dependent regimes 1,200 54–882 1–5 Rabosky ( 2014 ) and Mitchell and Rabosky ( 2016 )
Speciation rate evolves via diffusion process 1,200 25–1,208 1 Rabosky ( 2010 ), Beaulieu and O'Meara ( 2015 ), and Rabosky ( 2016 ) and this study

Size of diversification rate regimes might be an important factor in a tip rate metric's ability to accurately estimate rates. For example BAMM's statistical power in detecting a shift to a new rate regime is a function of the number of taxa in that rate regime, and tip rates for taxa from small regimes will more likely be parameterized according to the larger parent regime or the tree-wide average rate (Rabosky et al., 2017 ) this is the expected behaviour when BAMM fails to identify a rate shift. However, non-model-based approaches such as those examined in this study might be more accurate for small regimes. To explore how rate regime size influences the accuracy of tip rate metrics, we calculated the mean tip rate for each true rate regime from all multi-regime phylogenies (simulation datasets from Meyer & Wiens, 2017 Mitchell, Etienne, & Rabosky, 2018 Moore et al., 2016 Rabosky et al., 2017 ). We then calculated the Pearson correlation coefficient and the slope of a linear model between true and estimated mean regime rates. We explored the performance of all metrics with respect to regime sample size, as in Rabosky et al. ( 2017 , figure 13). For comparison, we repeated all performance summaries on tip rates estimated by applying a simple constant-rate birth–death (CRBD) process to each simulated phylogeny. This exercise is an important control, because it indicates how much error we would expect for each simulated phylogeny under the simplifying (incorrect) assumption that rates are constant among lineages and through time for each dataset.


This analysis is based on contributions to the RAINFOR network and database (, and supported by the Gordon and Betty Moore Foundation, National Environmental Research Council (grant numbers NE/I028122/1, NE/F005806/1), the European Commission [FP 5, 6 & 7 including the AMAZALERT (282664) and GEOCARBON (283080) projects)] and the Royal Society. Additional funding for fieldwork was provided by the National Geographic Society, Tropenbos International, European Commission, NASA Long-term Biosphere-Atmosphere Project in Amazonia (LBA), Brazilian National Council for Scientific and Technological Development (CNPq) including the long-term ecological research program CNPq/PELD Sítio 15 Transição Cerrado – Floresta Amazônica (558069/2009-6) and projeto INCT Processo 574008/2008-0, the National Institute for Amazonian Research (INPA), Brazil and the Tropical Ecology Assessment and Monitoring (TEAM) Network, a collaboration among Conservation International, the Missouri Botanical Garden, the Smithsonian Institution, and the Wildlife Conservation Society, and partially funded by these institutions, the Gordon and Betty Moore Foundation and Investissement d'Avenir grants of the French ANR (CEBA: ANR-10-LABX-0025 TULIP: ANR-10-LABX-0041), as well as other donors. OP, SL and GL-G are supported by the European Research Council project ‘Tropical Forests in the Changing Earth System’. We thank Tiina Saarkinen for assisting with clade ages for Fabaceae, David Greenberg for valuable discussions and Rick Condit for advice using lmerBayes.

Watch the video: Mitosis vs. Meiosis: Side by Side Comparison (May 2022).