Plant phylogenetic history explains in-stream decomposition at a global scale
Data Availability Statement: All raw and processed data, R-scripts, and figures are available at https://github.com/andrew-hipp/decomposition-phylogeny-2019.
Abstract
en- Evolutionary history and adaptation to climate shape plant traits. Some include leaf traits that influence litter quality. Thus, evolutionary history should affect litter decomposition, a crucial ecosystem process. In addition, litter decomposition is directly influenced by climate. We consequently expect plant phylogeny, adaptation and climate to jointly influence litter decomposition. These effects and their interactions have yet to be untangled at a global scale.
- Here we present an analysis of variation in litter decomposition rates in rivers and streams across 285 published studies for 239 species (from ferns to angiosperms) distributed at 494 locations world-wide. We estimated the relative contributions of climatic conditions and phylogenetic heritage on litter decomposition rates, partitioning phylogenetic from climatic effects at the site and species levels using phylogenetic eigenvector analysis and phylogenetic linear mixed models. In addition, we modelled transitions in decomposition rates under a suite of multiple adaptive-regime Ornstein–Uhlenbeck models to test the hypothesis that natural selection has shaped clade-level litter decomposition rates.
- Leaf litter decomposition rate exhibited a significant phylogenetic signal. Modelling decomposition rate as a function of location, climatic niche and phylogeny consistently recovered phylogeny alone as one of the top models in species-level analyses. Since many previous studies have focused on the same species across many locations, we also conducted analyses at the species × site level. Both phylogenetic and climatic factors were favoured in models of this analysis, but the single most important predictor for angiosperms and for all taxa analysed together was phylogeny alone.
- Synthesis. For plant species distributed globally at nearly 500 locations we found that plant phylogenetic history is a critically important predictor of litter decomposition rate in rivers and streams, explaining more of the variance in decomposition than site or climatic regime. Thus, our study demonstrates the influence of evolutionary history on suites of plant traits that shape a key ecosystem process.
Foreign Language Abstract Résumé
fr- L’histoire évolutive et l'adaptation au climat façonnent les caractéristiques des plantes. Le contrôle des traits foliaires sur la qualité des litières végétales pourrait sous-tendre l'existence de contraintes évolutives sur la décomposition des litières. Puisque ce processus écologique majeur est également directement modulé par le climat, il est important de mieux comprendre comment la phylogénie des plantes, leurs réponses adaptatives et le climat influencent ensemble la décomposition des litières à une échelle globale.
- Nous présentons ici les résultats d'une analyses quantitatives des données de décomposition de litières dans des rivières et des ruisseaux, issues de 285 publications couvrant 239 espèces de plantes (des fougères aux angiospermes) dans 494 lieux à travers le monde. Nous avons quantifié la contribution des conditions climatiques et de l'héritage phylogénétique aux variations du taux de décomposition des litières entre les espèces et entre les sites, à l'aide de régressions sur les vecteurs propres phylogénétiques et de modèles phylogénétiques linéaires à effets mixtes. Le modèle de Ornstein-Uhlenbeck a également été utilisé pour tenter de comprendre si et comment la sélection naturelle s'exerce sur la qualité des litières.
- Un signal phylogénétique cohérent a été détecté au sein des données analysées dans notre étude. La construction de modèles utilisant les informations géographiques, climatiques et phylogénétiques disponibles a permis d’établir la phylogénie comme un des meilleurs prédicteurs des variations interspécifiques du taux de décomposition des litières. Les variations spatiales de la décomposition étaient expliquées par les contraintes climatiques et phylogénétiques. Néanmoins, seule la phylogénie était nécessaire pour prédire les variations du taux de décomposition chez les angiospermes seulement ou bien sur l'ensemble des espèces de litière.
- Synthèses. En analysant les résultats d’études englobant un grand nombre d'espèces de plantes au sein de 500 sites environ, nous avons découvert que l'histoire évolutive des plantes prédisait mieux le taux de décomposition des litières dans les rivières et les ruisseaux que les données géographiques ou climatiques. Notre étude démontre comment les contraintes évolutives s'exerçant sur un ensemble de traits peuvent, de manière ultime, influencer un processus écologique majeur dans les écosystèmes.
1 INTRODUCTION
The past 20 years have seen a dramatic increase in our understanding of the influence of biodiversity on ecosystem functioning (Hooper et al., 2005; Cardinale et al., 2012; Gessner et al., 2010; Weisser et al., 2017). Expanding beyond the early emphasis of species richness effects on plant productivity, phylogenetic and genetic diversity have also been shown to influence ecosystem processes (Hughes, Inouye, Johnson, Underwood, & Vellend, 2008; Latzel et al., 2013; LeRoy, Whitham, Wooley, & Marks, 2007; Schweitzer et al., 2004; Sundqvist, Giesler, & Wardle, 2011). The exploration of the role of genetic patterns across species has generally taken one of two routes: (1) a diversity approach that uses phylogenetic diversity as a predictor of biodiversity–ecosystem function relationships (Flynn, Mirotchnick, Jain, Palmer, & Naeem, 2011; Srivastava, Cadotte, MacDonald, Marushia, & Mirotchnick, 2012), or (2) a phylogenetic approach that investigates how the evolution of traits and species relationships may influence ecosystem functioning (Cadotte et al., 2017; Donovan, Mason, Bowsher, Goolsby, & Ishibashi, 2014; Edwards, Still, & Donoghue, 2007). Studies taking the first approach have demonstrated that plant phylogenetic diversity influences plant productivity (Cadotte, Cardinale, & Oakley, 2008; Gravel et al., 2012), nutrient cycling (Cornelissen & Cornwell, 2014) and numerous species traits that subsequently affect ecosystem functioning (e.g. Cavender-Bares, Kozak, Fine, & Kembel, 2009; Díaz et al., 2013; Matthews et al., 2011; Senior et al., 2016). Despite a rapidly growing body of literature addressing these relationships, there have only been a few studies examining the influence of plant phylogenetic diversity on litter decomposition (Boyero et al., 2015; Cornwell et al., 2008; Liu et al., 2014; Makkonen et al., 2012; Pan et al., 2014). Even fewer studies have taken the second approach.
Only two previous studies in terrestrial ecosystems have explicitly investigated how plant phylogeny (and thus evolutionary history) shapes litter decomposition (Liu et al., 2014; Pan et al., 2014). Liu et al. (2014) applied a Brownian motion model to analyse data collected at a single site and found that mass loss was slower in basal angiosperms than in eudicot trees. Pan et al. (2014) compared three phylogenetic models for decomposition data also collected at a single site and found evidence that a constrained single-optimum Ornstein–Uhlenbeck (OU) model best fit their data, suggesting that decomposition rate is under stabilizing selection. No phylogenetic approaches have attempted to address the question of how important evolutionary history is in explaining litter decomposition rates across multiple sites. Understanding the relative contributions of contemporary environment and evolutionary heritage to decomposition rates globally has the potential to improve modelling carbon cycling on geological time-scales, by allowing us to estimate the decomposition rates characteristic of lineages that have variously dominated different habitats over deep time.
Plant litter decomposition is a crucial ecosystem process that drives carbon cycling in many terrestrial and aquatic environments. Shaded headwater streams make up the vast majority of river lengths globally (73.2%, Leopold, Wolman, & Miller, 1964) and strongly rely on plant litter inputs as a resource (Wallace, Eggert, Meyer, & Webster, 1997), relative to all lotic habitats. Plant litter decomposition at a global scale has been explained by both extrinsic (climate, latitude, altitude) and intrinsic (litter quality) factors (Aerts, 1997), but most previous studies have focused on terrestrial decomposition (Cornwell et al., 2008; Pietsch et al., 2014; Weedon et al., 2009). Results of these previous meta-analyses have generally shown that intrinsic species-level traits explain as much variation as extrinsic factors, if not more. Based on a review of 16 studies, Cornwell et al. (2008) found that species differences resulted in larger variation in decomposition rates than climatic differences. Specifically, litter traits like nitrogen and phosphorus concentrations led to faster decomposition and traits like lignin content, leaf mass area and the concentrations of water- and acid-soluble polysaccharides were associated with slower decomposition (Cornwell et al., 2008). Weedon et al. (2009) examined wood decomposition globally and found a phylogenetic distinction between angiosperms and gymnosperms with lower lignin, and higher N and P in the former, and a higher lignin-to-N ratio in the latter. Pietsch et al. (2014) explored the relationships between litter and wood traits and decomposition published in the two previous studies (Cornwell et al., 2008; Weedon et al., 2009) and found that wood and leaf decomposition rates were decoupled.
Despite decades of investigations into leaf litter decomposition in freshwaters, resulting in a wealth of information (Kaushik & Hynes, 1971; Petersen & Cummins, 1974; Tank, Rosi-Marshall, Griffiths, Entrekin, & Stephen, 2010; Webster & Benfield, 1986), previous global decomposition syntheses (Cornwell et al., 2008; Pietsch et al., 2014; Weedon et al., 2009) have until recently excluded aquatic studies. Recent global analyses of aquatic litter decomposition (Boyero et al., 2011, 2015; Follstad Shah et al., 2017; Handa et al., 2014; Tiegs et al., 2019) provide a strong framework for new comparative studies. The oldest of these (Boyero et al., 2011) compared breakdown rates at 26 globally distributed sites and found increased decomposition rates with increasing temperature and, when adjusted for temperature effects, with increasing latitude as a result of higher detritivore abundance. A follow-up study (Boyero et al., 2015) found an influence of phylogenetic diversity on breakdown rates of litter mixtures across 24 sites and 70 species. Follstad Shah et al. (2017) found that aquatic decomposition rates for hundreds of species across the globe tend to increase with latitude when adjusted for temperature effects using a large meta-analysis. Finally, using cotton strips as a standardized substrate at 514 stream and adjacent riparian sites world-wide, Tiegs et al. (2019) found wide variation in in-stream breakdown rates and decreased decomposition rates with increasing latitude.
However, none of these large-scale studies of aquatic leaf litter decomposition explored how the evolutionary history of individual species has shaped leaf litter decomposition. Phylogenetic diversity of plant communities has been addressed, but phylogenetic history of the constituent species integrates leaf trait evolution—indeed, the evolution of all plant traits—and as such may provide improved predictions of decomposition rates globally in linked terrestrial–aquatic ecosystems, where connections are mitigated by traits we may not predict a priori. Phylogenetic diversity of communities alone tells only part of the story of how evolutionary history shapes ecosystem function.
Here, we explore how both plant phylogenetic history and climatic variation among sites influence rates of leaf litter decomposition in streams distributed across the globe. Since climatic variables (mean, maximum and minimum temperature; precipitation; seasonality; isothermality) can influence both plant biogeography and the decomposition environment directly, we detangle these two influences (Dodds, Gido, Whiles, Daniels, & Grudzinski, 2015), testing the predictions that: (a) across broad spatial scales, decomposition rates will be strongly influenced by site × species interactions, with both factors influenced by climate, and (b) phylogenetic history will affect decomposition rates through its influences on species traits but also on the site × species interaction. As these predictions are not mutually exclusive, we tested them in a modelling framework that addresses the joint influences of phylogeny and climate at both levels of analysis, providing novel insights into how plant evolutionary history shapes stream ecosystem processes globally.
2 MATERIALS AND METHODS
2.1 Plant litter decomposition database compilation
Data compiled for this study were extracted from values reported in the scientific literature as detailed in Follstad Shah et al. (2017). Briefly, we used search terms “(leaf OR litter) AND (breakdown OR decomposition OR processing) AND (stream OR river)” to identify studies. This provided us with an initial list of articles published on May 13, 2011, which was updated by articles cited in review papers on leaf litter decomposition in streams by Follstad Shah et al., (2017). Our final list included 636 papers published between 1966 and 2011, of which 285 met our criteria that: (a) rates of leaf litter decomposition were measured in natural streams and rivers with perennial flow; and (b) stream temperature and decomposition rate constants were reported or could be calculated. We discarded studies including factors beyond local species differences (e.g. nutrient addition, predator exclusion, exotic species, long-distance reciprocal transplantations). However, we included in our dataset a variety of leaf litter decomposition methods used across studies, particularly litter bags of various size and mesh size, deployment periods and processing methods, to ensure a large sample size. Studies may have also included measurements of water chemistry, velocity, shredder abundance and microbial measures, but the inconsistency of these measurements made it impossible to comprehensively address these other factors.
For each article, we recorded the scientific name of the plant yielding the leaf litter and the corresponding decomposition rate constants as reported or calculated (see Follstad Shah et al., 2017 and methods therein), where missing, as 
, where mt is leaf litter mass remaining at time t and m0 is the initial leaf litter mass. Data compilation resulted in 3,189 records of leaf litter decomposition from 494 sites (Figure 1) for 239 plant taxa from 124 genera, 70 families and 34 orders. Decomposition rates ranged from 0.0002 to 0.7890 per day and resulted in a global, average in-stream decomposition rate of 0.0240 per day (Follstad Shah et al., 2017) with a 95% confidence interval of 0.0227–0.0253, and a median rate of 0.0127 per day.
Latitude and longitude data were extracted from each paper and mapped to identify potential errors which were subsequently rectified. Many of the studies lacked detailed data on environmental conditions (stream discharge, water chemistry, etc.), and a previous study exploring patterns in this dataset examined stream temperature (Follstad Shah et al., 2017), so we focused on broad climatic variables in this analysis. Climatic data were inferred for each site from the 10 min (~340 km2) resolution WorldClim v 2.0 dataset (Fick & Hijmans, 2017). All 19 of the WorldClim bioclimatic (“bioclim”) variables (https://www.worldclim.org/bioclim) were downloaded and used individually or summarized as ordination axes (see methods for Comparing environmental and phylogenetic drivers below). R scripts for executing analyses and downloading data have been archived at (https://github.com/andrew-hipp/decomposition-phylogeny-2019).
2.2 Plant phylogenetic data organization
We used Phylomatic (Beaulieu, Ree, Cavender-Bares, Weiblen, & Donoghue, 2012; Webb & Donoghue, 2005) to assemble a base phylogenetic tree for all species for which we had litter decomposition rate data (Appendix S1), using Phylomatic base tree R20120829 as the starting megatree, based on Angiosperm Phylogeny Group III (Haston, Richardson, Stevens, Chase, & Harris, 2009; Qian & Jin, 2016). We normalized scientific names using the Taxonomic Name Resolution Service (Boyle et al., 2013) and corrected as needed with names used in published phylogenies (https://github.com/andrew-hipp/decomposition-phylogeny-2019). Mesquite (Maddison & Maddison, 2018) was used to manually resolve branches within families, as the supertree is not fully resolved for all families, based on those same phylogenetic studies (Appendix S1). Node ages of the resulting tree were then calibrated using more recent angiosperm phylogenies (Bell, Soltis, & Soltis, 2010) and the simple branch length adjuster tool (BLADJ) to even out node spacing between calibration points. Initial analyses conducted on the tree with unmodified ages did not suggest different interpretations than analyses conducted on the Phylomatic tree with updated ages, suggesting that our results are robust to a range of branch length assumptions. Throughout the paper, we report analyses conducted on the tree with revised clade ages and results for the entire dataset as well as a dataset including just the angiosperms. The tree was visualized and exported for publication in R using the ggplot2 and ggtree packages (Yu et al., 2016) and custom scripts (https://github.com/andrew-hipp/decomposition-phylogeny-2019). Site-level phylogenetic diversity was calculated for all sites where more than one taxon occurred using mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) in picante v. 1.7 for R (Kembel et al., 2010).
For analyses at the site level using generalized linear mixed models, we created a composite tree with one tip per species per site, creating a polytomy for each species with a depth of 1.5 Ma to ensure that all species polytomies were younger than the youngest node in the phylogeny. Thus, each species × site combination had one tip in the tree, allowing us to partition the site and phylogenetic contributions to variance in litter decomposition rates.
2.3 Phylogenetic signal and phylogenetic transitions on decomposition rate
We estimated phylogenetic signal, defined as the “tendency for related species to resemble each other more than they resemble species drawn at random from the tree”, (Blomberg & Garland, 2002) in litter decomposition rates using Blomberg's K (Blomberg, Garland, & Ives, 2003). Blomberg's K typically scales from 0 for a character with no phylogenetic autocorrelation to 1 for a character that has phylogenetic autocorrelation (i.e. similarity in trait value between close relatives) equivalent to expectations for a trait evolving according to a Brownian motion (random walk) process on the tree being evaluated (Figure 2). K may range to greater than 1 for traits that exhibit greater phylogenetic autocorrelation than expected under the assumption of Brownian motion. Stated another way, K > 1 indicates that close relatives are more similar than expected if the trait being studied evolves according to a random walk. Significance was assessed by permuting tip states 4,999 times, calculating K for each permutation, then calculating p as 2 × the rank position of Kobserved; thus yielding a two-tailed p-value, reflecting the fact that both clustering (relatives more similar than expected) and overdispersion (close relatives less similar than expected) are possible evolutionary outcomes. These analyses were conducted with the package picante 1.7 (Kembel et al., 2010) in r 3.4.4 (R Core Team, 2018). For comparison with previous studies, we also calculated Pagel's lambda (Münkemüller et al., 2012; Pagel, 1999), a scalar of the off-diagonal cells of the covariance matrix estimated from the tree, which scales from 0 to 1 (or a bit higher, depending on the tree structure). At λ = 0, the covariance is equivalent to a star-shaped phylogeny, where all species are equally related, reflecting the case where phylogeny has no effect on litter decomposition rates. At λ = 1, the covariance between any pair of tips is predicted directly by the relative distance from the root of the tree to the most recent common ancestor of those tips; this models a Brownian motion process in which the branch lengths of the tree are known without error. The λ model was fitted using the fitContinuous function in the package geiger 2.0.6 (Harmon, Weir, Brock, Glor, & Challenger, 2008) in r 3.4.4.
We identified phylogenetic shifts in litter decomposition rate by evaluating the relative support for alternative OU models, which model transitions in trait values as responses to shifting selective regimes (Butler & King, 2004; Hansen, 1997; Martins & Hansen, 1997). Analysis was performed using the l1ou + IC method (Khabbazian, Kriebel, Rohe, & Ané, 2016), a stepwise model-selection and evaluation approach implemented using the least absolute shrinkage and selection operator (lasso) (Tibshirani, 1996, 2011) method and a modified information criterion that accounts for clade sizes entailed by shifts in litter decomposition rate, modelled as shifting optima in trait space (where the trait is litter decomposition rate). We compared results with an Expectation Maximization (EM) search algorithm (Bastide, Ané, Robin, & Mariadassou, 2018), which generally searches parameter space more rapidly. Analyses were conducted in the l1ou (Khabbazian et al., 2016) and PhylogeneticEM (Bastide, Mariadassou, & Robin, 2017) packages of R 3.4.4.
We reconstructed phylogenetic transitions in the rate of evolution of litter decomposition rate (i.e. the rate of evolution of log(k)) using reversible-jump Markov-chain Monte Carlo (rjMCMC; Green, 1995) as implemented in the rjmcmc.bm function of geiger 2.0.6 (Eastman, Alfaro, Joyce, Hipp, & Harmon, 2011; Harmon et al., 2008). In this method, every branch on the phylogeny was assigned to a class of rates, corresponding to presumed shifts in the rate of evolution on the phylogenetic tree. The relative rates of all branches within a taxonomic class were allowed to vary. The assignments of branches to classes were also allowed to vary by allowing transitions in rate to appear or disappear from the tree.
2.4 Comparing environmental and phylogenetic drivers
We used two methods to compare phylogenetic and climatic models to predict global patterns of litter decomposition rates. First, we used an approach based on phylogenetic eigenvector regression (PVR; Diniz-Filho, Sant’Ana, & Bini, 1998) to partition variance in litter decomposition rate into (a) an environmental component based on the four strongest bioclim variables (Bioclim1: annual mean temperature; Bioclim4: temperature seasonality; Bioclim12: annual precipitation; and Bioclim14: precipitation of the driest month); (b) four non-metric multidimensional scaling (NMDS) ordination axes created using a Euclidean distance measure from all 19 bioclim variables (https://www.worldclim.org/bioclim) (and for comparison, four principle component analysis (PCA) axes created using Euclidean distance measures with all 19 variables); (c) a phylogenetic component based on the first four eigenvectors from principal coordinate analysis (PCoA) of a pairwise phylogenetic distance matrix; and (d) site identity for analyses conducted at the site level. Site identity was added as a random effect to account for any aspects of site that are not directly measured; in this way, site identity in our model is analogous to plot or block in a designed experiment. A stepwise model selection process is generally employed to select phylogenetic eigenvectors for analysis (Diniz-Filho et al., 2012; Gonçlaves-Souza, Diniz-Filho, & Romero, 2014). However, the lack of an explicit evolutionary model underlying phylogenetic eigenvector analysis (Freckleton, Cooper, & Jetz, 2011) combined with the relative arbitrariness of eigenvector selection—selecting all 238 eigenvectors would overfit the data, while selecting anything less than all fails to capture all phylogenetic data—makes PVR an approximate approach. Consequently, we used only four phylogenetic eigenvectors, which cumulatively explain 55% of the total variance in the phylogeny when the relative eigenvalues are used as estimates of the variance explained by each PCoA axis. Axes 2 and 3 of the angiosperm ordination recovered differences among broadly overlapping clades; consequently, we used PCoA axes 1, 2, 4 and 5 for model evaluation in the angiosperms by species models, but this had a negligible effect on model ranking and support. These four eigenvectors were included in models with either four bioclim variables or four bioclim NMDS axes for a comparison of the relative predictive value of both climate and phylogeny. We selected the same number of axes, four, from both the phylogeny PCoA and the climatic NMDS to avoid biasing our linear regressions towards detecting either climate or phylogeny as a better predictor of decomposition rate. In both cases, the axes explain a fair amount of variation in the data. For the phylogenetic PCoA, the first four axes explain 33.0, 8.6, 8.0 and 5.4% of variation respectively. For NMDS, variance partitioning is not the target of each axis, but ordination stress—a measure of the mismatch between the K-dimensional ordination and the distance matrix used to generate it. Ordination stress drops from 0.214 at K = 1 to 0.035 at K = 4; all dimensions above K = 4 reduce stress by less than 0.02 per increase in dimension. To assess the effect of ordination method on model support and parameter estimates, we replicated the analyses with the first four axes of a PCA using the 19 bioclim climatic predictors and found that there is no difference in model support and negligible difference in the estimated regression parameters (Supplemental Table S1). Code to do all analyses is in the repository for this paper (https://github.com/andrew-hipp/decomposition-phylogeny-2019). Linear regression models were compared based on Akaike's Information Criterion (AIC; Akaike, 1974) and coefficients of determination (R2), and model confidence intervals were constructed using cumulative AIC weights (Burnham & Anderson, 2002; Burnham, Anderson, & Huyvaert, 2011). All analyses were conducted using the ape (Paradis, Claude, & Strimmer, 2004) and vegan (Oksanen, Blanchet, Kindt, Legendre, & O’Hara, 2016) packages of R and custom scripts available at https://github.com/andrew-hipp/decomposition-phylogeny-2019.
Our second approach used phylogenetic generalized linear mixed models (PGLMM; Hadfield & Nakagawa, 2010) in a Bayesian framework to quantify the effects of phylogeny, climate, and site on litter decomposition rate more precisely. This approach uses the phylogenetic covariance matrix as a predictor of covariance in trait value among tips, and it scales the importance of the phylogeny by rescaling the covariance matrix. For a species-level analysis, this is analogous to fitting a phylogenetic generalized least squares model while simultaneously estimating a scalar of the off-diagonal elements of the covariance matrix (i.e. Pagel's λ) (Housworth, Martins, & Lynch, 2004; Pagel, 1997, 1999; Revell, 2010). In species-level analyses, we treat phylogeny as a random effect and bioclim variables as fixed effects. For a site-level analysis, PGLMM offers the flexibility of treating both phylogeny and site as random effects and bioclim variables as fixed effects. We compared PGLMM models based on Deviance Information Criterion (DIC; van der Linde, 2005; Spiegelhalter, Best, Carlin, & Linde, 2002) and residual model variance, and model intervals were constructed using DIC weights. Analyses were conducted in R using the MCMCglmm package (Hadfield, 2010) and custom scripts available at https://github.com/andrew-hipp/decomposition-phylogeny-2019.
3 RESULTS
Litter decomposition rates (k/day) in stream environments showed a wide range of variation across plant families (min–max: 0.0036–0.1522) and plant orders globally (Table 1). Our results (averaged by plant family) are compared to an early synthesis of decomposition rates (Webster & Benfield, 1986; Figure 3) to provide a more complete picture of global decomposition across plant families.
| Clade | Order | Family | N | mean k/day | SD | SEM |
|---|---|---|---|---|---|---|
| Fern/fern allies | eupolypod II | Blechnaceae | 1 | 0.0066 | – | – |
| Fern/fern allies | Marattiales | Marattiaceae | 1 | 0.013 | – | – |
| Gymnosperm | Cupressales | Cupressaceae | 1 | 0.008 | – | – |
| Gymnosperm | Cupressales | Podocarpaceae | 1 | 0.016 | – | – |
| Gymnosperm | Pinales | Pinaceae | 8 | 0.0058 | 0.0054 | 0.0019 |
| Monocot | Alismatales | Potamogetonaceae | 2 | 0.0091 | 0.0077 | 0.0054 |
| Monocot | Arecales | Arecaceae | 1 | 0.008 | – | – |
| Monocot | Commelinales | Pontederiaceae | 1 | 0.0122 | – | – |
| Monocot | Pandanales | Pandanaceae | 1 | 0.0144 | – | – |
| Monocot | Poales | Cyperaceae | 3 | 0.0127 | 0.011 | 0.0063 |
| Monocot | Poales | Poaceae | 12 | 0.0108 | 0.0086 | 0.0025 |
| Monocot | Poales | Typhaceae | 1 | 0.0252 | – | – |
| Magnoliid | Laurales | Atherospermataceae | 1 | 0.0072 | – | – |
| Magnoliid | Laurales | Lauraceae | 4 | 0.006 | 0.0043 | 0.0022 |
| Magnoliid | Magnoliales | Magnoliaceae | 1 | 0.0155 | – | – |
| Eudicot | Apiales | Araliaceae | 2 | 0.0291 | 0.01 | 0.0071 |
| Eudicot | Apiales | Pittosporaceae | 2 | 0.0296 | 0.0195 | 0.0138 |
| Eudicot | Aquifoliales | Aquifoliaceae | 2 | 0.0303 | 0.0293 | 0.0207 |
| Eudicot | Asterales | Asteraceae | 1 | 0.0216 | – | – |
| Eudicot | Brassicales | Brassicaceae | 1 | 0.0845 | – | – |
| Eudicot | Caryophyllales | Polygonaceae | 2 | 0.0078 | 0.0032 | 0.0023 |
| Eudicot | Caryophyllales | Tamaricaceae | 2 | 0.0064 | 0.0022 | 0.0016 |
| Eudicot | Cornales | Cornaceae | 3 | 0.04 | 0.0274 | 0.0158 |
| Eudicot | Cornales | Nyssaceae | 1 | 0.0313 | – | – |
| Eudicot | Dipsacales | Caprifoliaceae | 1 | 0.1018 | – | – |
| Eudicot | Ericales | Ericaceae | 3 | 0.004 | 0.0032 | 0.0019 |
| Eudicot | Ericales | Primulaceae | 1 | 0.007 | – | – |
| Eudicot | Fabales | Fabaceae | 6 | 0.0311 | 0.035 | 0.0143 |
| Eudicot | Fagales | Betulaceae | 17 | 0.0229 | 0.0233 | 0.0057 |
| Eudicot | Fagales | Casuarinaceae | 1 | 0.0209 | – | – |
| Eudicot | Fagales | Fagaceae | 18 | 0.0114 | 0.0055 | 0.0013 |
| Eudicot | Fagales | Juglandaceae | 3 | 0.0095 | 0.0039 | 0.0022 |
| Eudicot | Fagales | Myricaceae | 1 | 0.039 | – | – |
| Eudicot | Fagales | Nothofagaceae | 4 | 0.0106 | 0.0071 | 0.0035 |
| Eudicot | Gentianales | Apocynaceae | 2 | 0.0237 | 0.0181 | 0.0128 |
| Eudicot | Gentianales | Rubiaceae | 2 | 0.0476 | 6.00E-04 | 4.00E-04 |
| Eudicot | Icacinales | Icacinaceae | 1 | 0.006 | – | – |
| Eudicot | Lamiales | Oleaceae | 4 | 0.0182 | 0.0085 | 0.0042 |
| Eudicot | Lamiales | Scrophulariaceae | 1 | 0.1522 | – | – |
| Eudicot | Lamiales | Verbenaceae | 1 | 0.0182 | – | – |
| Eudicot | Malpighiales | Euphorbiaceae | 5 | 0.0411 | 0.023 | 0.0103 |
| Eudicot | Malpighiales | Salicaceae | 28 | 0.0158 | 0.0133 | 0.0025 |
| Eudicot | Malpighiales | Violaceae | 1 | 0.0477 | – | – |
| Eudicot | Malvales | Malvaceae | 5 | 0.0858 | 0.1407 | 0.0629 |
| Eudicot | Myrtales | Combretaceae | 2 | 0.025 | 0.0198 | 0.014 |
| Eudicot | Myrtales | Lythraceae | 1 | 0.0333 | – | – |
| Eudicot | Myrtales | Melastomataceae | 2 | 0.0064 | 0.0027 | 0.0019 |
| Eudicot | Myrtales | Myrtaceae | 20 | 0.0115 | 0.0091 | 0.002 |
| Eudicot | Myrtales | Onagraceae | 1 | 0.0165 | – | – |
| Eudicot | Oxalidales | Cunoniaceae | 2 | 0.0363 | 0.0471 | 0.0333 |
| Eudicot | Oxalidales | Elaeocarpaceae | 2 | 0.0286 | 0.0369 | 0.0261 |
| Eudicot | Proteales | Platanaceae | 6 | 0.0057 | 0.0019 | 8.00E-04 |
| Eudicot | Proteales | Proteaceae | 2 | 0.009 | 0.0022 | 0.0016 |
| Eudicot | Ranunculales | Eupteleaceae | 1 | 0.0133 | – | – |
| Eudicot | Rosales | Elaeagnaceae | 1 | 0.024 | – | – |
| Eudicot | Rosales | Moraceae | 6 | 0.0388 | 0.0308 | 0.0126 |
| Eudicot | Rosales | Rhamnaceae | 1 | 0.0286 | – | – |
| Eudicot | Rosales | Rosaceae | 4 | 0.0394 | 0.0443 | 0.0222 |
| Eudicot | Rosales | Ulmaceae | 7 | 0.0705 | 0.0727 | 0.0275 |
| Eudicot | Rosales | Urticaceae | 1 | 0.02 | – | – |
| Eudicot | Sapindales | Anacardiaceae | 2 | 0.0202 | 0.0078 | 0.0056 |
| Eudicot | Sapindales | Burseraceae | 2 | 0.0036 | 4.00E-04 | 2.00E-04 |
| Eudicot | Sapindales | Hippocastanaceae | 1 | 0.0066 | – | – |
| Eudicot | Sapindales | Meliaceae | 1 | 0.023 | – | – |
| Eudicot | Sapindales | Rutaceae | 1 | 0.019 | – | – |
| Eudicot | Sapindales | Sapindaceae | 8 | 0.0232 | 0.0146 | 0.0052 |
| Eudicot | Sapindales | Simaroubaceae | 1 | 0.0945 | – | – |
| Eudicot | Saxifragales | Haloragaceae | 1 | 0.025 | – | – |
| Eudicot | Saxifragales | Hamamelidaceae | 2 | 0.0172 | 0.0029 | 0.002 |
| Eudicot | Vitales | Vitaceae | 1 | 0.0142 | – | – |
Note
- Standard deviation (SD) and standard error (SEM) are provided for families with N > 1.
3.1 Phylogenetic signal and phylogenetic transitions in litter decomposition rate
The litter decomposition rate data exhibited a significant phylogenetic signal in both the all-taxa tree (Figure 4; K = 0.177, p = .01; λ = 0.977) and the angiosperm tree (K = 0.376, p = .004; λ = 0.961; Table 2). Because plants are distributed non-randomly across the globe, climatic variables averaged for each species all exhibited some phylogenetic signal. Mean annual temperature (BIO1) and temperature seasonality (BIO4) each exhibited significant phylogenetic signal (on the angiosperm tree, K = 0.230–0.266, p ≤ .002, λ = 0.59–0.70; on the all-taxa tree, K = 0.097–0.117, p ≤ .006, λ = 0.74–0.83), and precipitation of the driest month (BIO14) exhibited significant phylogenetic signal on the angiosperm tree (K = 0.227, p = .038, λ = 0.23), but a non-significant phylogenetic signal on the all-taxa tree (K = 0.1, p = .054, λ = 0.0), while mean annual precipitation (BIO12) exhibited weak and, in our dataset, non-significant phylogenetic signal (K = 0.094–0.208, p = .078–0.112, λ = 0.51–0.64; Table 2). Ordination of these data in four dimensions (Figure S1) was well-supported (stress = 0.035), and axes 1 and 3 exhibited weak but significant phylogenetic signal on the angiosperm tree (K = 0.200–0.255, p < .042, λ = 0.35–0.70) and the all-taxa tree (K = 0.093–0.111, p < .032, λ = 0.52–0.83).
| Angiosperm tree | All taxa tree | |
|---|---|---|
| k (day−1) | K = 0.376, p = 0.004, λ = 0.961 | K = 0.177, p = 0.010, λ = 0.977 |
| ln k (day−1) | K = 0.226, p = 0.010, λ = 0.536 | K = 0.122, p = 0.004, λ = 0.716 |
| Bio1 | K = 0.230, p = 0.002, λ = 0.587 | K = 0.097, p = 0.006, λ = 0.743 |
| Bio12 | K = 0.208, p = 0.078, λ = 0.505 | K = 0.094, p = 0.122, λ = 0.638 |
| Bio4 | K = 0.266, p < 0.001, λ = 0.703 | K = 0.117, p < 0.001, λ = 0.833 |
| Bio14 | K = 0.227, p = 0.038, λ = 0.225 | K = 0.100, p = 0.054, λ = 0.000 |
| NMDS1 | K = 0.255, p = 0.002, λ = 0.697 | K = 0.111, p < 0.001, λ = 0.826 |
| NMDS2 | K = 0.186, p = 0.204, λ = 0.106 | K = 0.080, p = 0.262, λ = 0.000 |
| NMDS3 | K = 0.200, p = 0.042, λ = 0.343 | K = 0.093, p = 0.032, λ = 0.510 |
| NMDS4 | K = 0.192, p = 0.058, λ = 0.000 | K = 0.090, p = 0.072, λ = 0.000 |
Analyses of alternative multiple-regime OU models identified a significant increase in litter decomposition rate at the base of the eudicots (Figure S2). Both inferred a relatively high rate of adaptation, corresponding to 5.0%–5.5% of the total tree length (for the phylogenetic lasso, α = 14.36 [t½ = 0.0483] on a tree arbitrarily rescaled to height 1.0; for the E-M algorithm, α = 0.0305 [t½ = 22.73] on a tree spanning 432 million years, t½ = 5.26% of total tree depth). While interpreting the rate of adaptation (α) is not perfectly straightforward (Bastide et al., 2018), this rate of adaptation suggests that the transition to a higher litter decomposition rate at the base of the eudicots occurred too rapidly to be explained by a purely random walk model, that is, by the variance we expected as lineages diverge from one another (Felsenstein, 1985). This in turn suggests that natural selection on leaf traits shaped this transition, driving accelerated rates of litter decomposition in the angiosperms. By contrast, analysis of evolutionary rates using rjMCMC suggests that the evolution of litter decomposition rates is relatively homogeneous: no lineages were recovered as evolving at a significantly higher or lower rate than other clades (Figure S3).
3.2 Comparing environmental and phylogenetic drivers
Top-ranked linear models for both the species-level and site-level analyses included phylogenetic PCoA axes, irrespective of whether the all-taxa phylogeny or angiosperms-only phylogeny was analysed (Table 3). Climatic predictors, whether included directly in the model as bioclim variables, or indirectly as the first four NMDS axes, appeared in the 95% model confidence set (based on cumulative AIC weight) only in combination with phylogeny, whereas phylogeny alone had the lowest AIC of any of the individual models in the species-level analysis of the angiosperm dataset. The first axis of the phylogenetic PCoA had a stronger effect than climatic predictors in multiple regression as estimated using regression coefficients on rescaled data for the all-taxa analyses and the angiosperm dataset analysed by species. For the angiosperm dataset analysed by site, climatic predictors had a stronger effect (Table 3). Among all models in the 95% confidence set, mean temperature (BIO1) had the strongest climatic effect on litter decomposition rate as estimated by normalized partial regression coefficients, followed by either mean annual precipitation (BIO12) or temperature seasonality (BIO4) (Table 3).
| Taxon set | Analysis level | Model | R 2 | AIC | delta AIC | Sum AICw | bio1 | bio12 | bio14 | bioclim4 | NMDS1 | NMDS2 | NMDS3 | NMDS4 | phylo Axis1 | phylo Axis2 | phylo Axis3 | phylo Axis4 | phylo Axis5 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| allTaxa | bySp | NMDS & phyloPCOA | 0.114 | 668.205 | 0 | 0.394 | – | – | – | – | 0.137 (p = .033) | 0.006 (p = .929) | 0.075 (p = .238) | −0.105 (p = .093) | −0.213 (p = .001) | 0.13 (p = .042) | 0.113 (p = .074) | 0.022 (p = .727) | – |
| allTaxa | bySp | bioclim & phyloPCOA | 0.114 | 668.302 | 0.097 | 0.77 | 0.129 (p = .224) | 0.164 (p = .205) | −0.008 (p = .939) | 0.122 (p = .198) | – | – | – | – | −0.22 (p = .001) | 0.13 (p = .043) | 0.118 (p = .062) | 0.019 (p = .771) | – |
| allTaxa | bySp | phyloPCOA | 0.08 | 669.303 | 1.098 | 0.998 | – | – | – | – | – | – | – | – | −0.236 (p < .001) | 0.112 (p = .075) | 0.101 (p = .109) | 0.038 (p = .544) | – |
| allTaxa | bySp | NMDS | 0.04 | 679.394 | 11.189 | 0.999 | – | – | – | – | 0.143 (p = .026) | −0.005 (p = .937) | 0.051 (p = .428) | −0.131 (p = .041) | – | – | – | – | – |
| allTaxa | bySp | bioclim | 0.035 | 680.678 | 12.473 | 1 | 0.13 (p = .222) | 0.131 (p = .325) | 0.032 (p = .775) | 0.11 (p = .256) | – | – | – | – | – | – | – | – | – |
| allTaxa | bySpxSite | bioclim & phyloPCOA | 0.096 | 2,523.736 | 0 | 1 | 0.131 (p = .012) | 0.047 (p = .41) | 0.002 (p = .958) | −0.077 (p = .111) | – | – | – | – | −0.197 (p < .001) | 0.032 (p = .39) | 0.046 (p = .206) | −0.039 (p = .23) | – |
| allTaxa | bySpxSite | bioclim | 0.048 | 2,562.93 | 39.194 | 1 | 0.12 (p = .019) | 0.05 (p = .387) | −0.015 (p = .757) | −0.084 (p = .085) | – | – | – | – | – | – | – | – | – |
| allTaxa | bySpxSite | phyloPCOA | 0.047 | 2,563.876 | 40.14 | 1 | – | – | – | – | – | – | – | – | −0.202 (p < .001) | 0.008 (p = .839) | 0.045 (p = .225) | 0.001 (p = .978) | – |
| allTaxa | bySpxSite | NMDS & phyloPCOA | 0.052 | 2,566.611 | 42.875 | 1 | – | – | – | – | 0.074 (p = .031) | −0.029 (p = .395) | 0.011 (p = .733) | −0.014 (p = .664) | −0.203 (p < .001) | 0.012 (p = .748) | 0.048 (p = .203) | −0.015 (p = .66) | – |
| allTaxa | bySpxSite | NMDS | 0.005 | 2,603.493 | 79.757 | 1 | – | – | – | – | 0.068 (p = .044) | −0.019 (p = .578) | −0.001 (p = .984) | −0.016 (p = .635) | – | – | – | – | – |
| angiosperms | bySp | NMDS & phyloPCOA | 0.082 | 643.89 | 0 | 0.356 | – | – | – | – | 0.149 (p = .027) | 0.003 (p = .963) | 0.055 (p = .421) | −0.103 (p = .116) | −0.249 (p = .088) | −0.03 (p = .837) | – | 0.023 (p = .729) | −0.04 (p = .555) |
| angiosperms | bySp | phyloPCOA | 0.047 | 644.313 | 0.423 | 0.643 | – | – | – | – | – | – | – | – | −0.249 (p = .09) | −0.047 (p = .747) | – | 0.044 (p = .506) | −0.035 (p = .591) |
| angiosperms | bySp | bioclim & phyloPCOA | 0.08 | 644.384 | 0.494 | 0.921 | 0.13 (p = .245) | 0.122 (p = .387) | 0.023 (p = .842) | 0.087 (p = .396) | – | – | – | – | −0.252 (p = .085) | −0.03 (p = .837) | – | 0.018 (p = .786) | −0.033 (p = .623) |
| angiosperms | bySp | NMDS | 0.031 | 648.027 | 4.137 | 0.966 | – | – | – | – | 0.13 (p = .051) | −0.031 (p = .639) | 0.04 (p = .545) | −0.108 (p = .102) | – | – | – | – | – |
| angiosperms | bySp | bioclim | 0.029 | 648.583 | 4.693 | 1 | 0.074 (p = .505) | 0.132 (p = .341) | 0.036 (p = .757) | 0.072 (p = .476) | – | – | – | – | – | – | – | – | – |
| angiosperms | bySpxSite | bioclim & phyloPCOA | 0.076 | 2,484.951 | 0 | 1 | 0.156 (p = .004) | 0.032 (p = .576) | 0.024 (p = .618) | −0.076 (p = .129) | – | – | – | – | 0.068 (p = .045) | −0.061 (p = .071) | 0.123 (p < .001) | 0.085 (p = .009) | – |
| angiosperms | bySpxSite | bioclim | 0.047 | 2,504.595 | 19.644 | 1 | 0.089 (p = .089) | 0.04 (p = .499) | 0.002 (p = .959) | −0.117 (p = .018) | – | – | – | – | – | – | – | – | – |
| angiosperms | bySpxSite | phyloPCOA | 0.02 | 2,529.733 | 44.782 | 1 | – | – | – | – | – | – | – | – | 0.046 (p = .175) | −0.011 (p = .745) | 0.087 (p = .011) | 0.095 (p = .004) | – |
| angiosperms | bySpxSite | NMDS & phyloPCOA | 0.027 | 2,531.992 | 47.041 | 1 | – | – | – | – | 0.083 (p = .022) | −0.022 (p = .523) | 0.002 (p = .945) | −0.031 (p = .382) | 0.054 (p = .132) | −0.032 (p = .359) | 0.103 (p = .004) | 0.087 (p = .01) | – |
| angiosperms | bySpxSite | NMDS | 0.005 | 2,543.133 | 58.182 | 1 | – | – | – | – | 0.064 (p = .061) | −0.05 (p = .145) | 0.0 (p = .994) | 0.001 (p = .975) | – | – | – | – | – |
Note
- All predictors used in the models are indicated by the presence of a regression coefficient; parameters not included in each model are indicated by blanks in the table. Models were run using two sets of data: (a) all taxa in the dataset, and (b) just the angiosperms in the dataset. Analyses were conducted by species and by species × site. Results include coefficients of determination (R2), Akaike's Information Criterion (AIC), differences between the top-ranked model and others (delta AIC), and the cumulative AIC weight (Sum AICw), with lowest AIC scores reported at the top. Regression coefficients are all estimated on data rescaled to unit variance and mean of 0 and are highlighted in bold when p-values are less than 0.05. "NMDS" in the model name indicates that nonmetric multidimensional scaling axes for the Bioclim data were used in lieu of individual Bioclim variables.
Including phylogeny as a random effect in analyses by species reduced the residual variance by approximately 40% in the linear mixed models (Table 4). Including just site as a random effect reduced residual variance from 1.0 to 0.58–0.59 in analyses conducted by site; including just phylogeny (also as a random effect) reduced residual variance from 1.0 to 0.62 or 0.65 (Table 4). The single most important predictor in both sets of taxa and both analysis levels was phylogeny, which was for all sets of analyses the lowest-DIC single-predictor model (Table 4). Including climate as a suite of fixed effects did not improve model fit over that of the model including only phylogeny as a random effect in analyses conducted by species. In analyses conducted by site, the model including site and phylogeny as random effects and climate as a suite of fixed effects was by far the best fit model, with phylogeny explaining ca. 6 × the variance explained by site in the all taxa dataset, and ca. 2 × the variance explained by site in the angiosperm dataset (Table 4, site var. and phylo var. columns). In both datasets, mean annual temperature (BIO1) had approximately twice the effect size of mean annual rainfall (BIO12), but neither effect was significant (Table 4). Using ordination to determine the influence of climatic factors on species-specific decomposition rates (239 species plotted in climatic space, Figure 5), demonstrated that climate is a relatively poor predictor of decomposition rates at a global scale.
| Taxon set | Analysis level | Model | DIC | delta DIC | Sum DICw | Residual var | Site var | Phylo var | bio1 | bio12 | bio4 | bio14 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| allTaxa | bySp | phyloPCOA | 587.107 | 0 | 0.977 | 0.524 | – | 1.42 | – | – | – | – |
| allTaxa | bySp | NMDS & phyloPCOA | 594.644 | 7.537 | 1 | 0.542 | – | 1.31 | 0.002 (p = .238) | 0.000 (p = .972) | 0.000 (p = .471) | 0.001 (p = .591) |
| allTaxa | bySp | NMDS | 659.939 | 72.832 | 1 | 0.907 | – | – | 0.002 (p = .228) | 0.000 (p = .323) | 0.000 (p = .258) | 0.001 (p = .774) |
| allTaxa | bySp | no Predictors | 660.457 | 73.35 | 1 | 0.922 | – | – | – | – | – | – |
| allTaxa | bySpxSite | phyloPCOA & NMDS & site | 1612 | 0 | 1 | 0.203 | 0.394 | 2.263 | 0.122 (p = .045) | −0.053 (p = .436) | −0.125 (p = .050) | −0.014 (p = .822) |
| allTaxa | bySpxSite | phyloPCOA & site | 1627.434 | 15.434 | 1 | 0.206 | 0.414 | 2.182 | – | – | – | – |
| allTaxa | bySpxSite | phyloPCOA & NMDS | 2,193.091 | 581.091 | 1 | 0.563 | – | 2.37 | 0.182 (p = .000) | −0.016 (p = .786) | −0.121 (p = .026) | 0.008 (p = .868) |
| allTaxa | bySpxSite | phyloPCOA | 2,274.177 | 662.177 | 1 | 0.627 | – | 1.92 | – | – | – | – |
| allTaxa | bySpxSite | NMDS & site | 2,341.481 | 729.481 | 1 | 0.587 | 0.354 | – | 0.121 (p = .036) | −0.009 (p = .905) | −0.109 (p = .078) | −0.012 (p = .832) |
| allTaxa | bySpxSite | site | 2,354.479 | 742.479 | 1 | 0.591 | 0.38 | – | – | – | – | – |
| allTaxa | bySpxSite | NMDS | 2,562.968 | 950.968 | 1 | 0.956 | – | – | 0.120 (p = .018) | 0.050 (p = .409) | −0.083 (p = .084) | −0.013 (p = .792) |
| allTaxa | bySpxSite | no Predictors | 2,599.644 | 987.644 | 1 | 1.002 | – | – | – | – | – | – |
| Angiosperms | bySp | phyloPCOA | 574.868 | 0 | 0.939 | 0.559 | – | 0.681 | – | – | – | – |
| Angiosperms | bySp | NMDS & phyloPCOA | 580.337 | 5.469 | 1 | 0.569 | – | 0.68 | 0.102 (p = .324) | −0.022 (p = .854) | 0.038 (p = .683) | 0.068 (p = .516) |
| Angiosperms | bySp | no Predictors | 640.431 | 65.563 | 1 | 0.978 | – | – | – | – | – | – |
| Angiosperms | bySp | NMDS | 641.865 | 66.997 | 1 | 0.966 | – | – | 0.073 (p = .488) | 0.129 (p = .346) | 0.071 (p = .462) | 0.035 (p = .780) |
| Angiosperms | bySpxSite | phyloPCOA & NMDS & site | 1604.314 | 0 | 0.999 | 0.209 | 0.41 | 0.98 | 0.111 (p = .075) | −0.050 (p = .478) | −0.139 (p = .041) | −0.019 (p = .739) |
| Angiosperms | bySpxSite | phyloPCOA & site | 1619.028 | 14.714 | 1 | 0.213 | 0.431 | 0.934 | – | – | – | – |
| Angiosperms | bySpxSite | phyloPCOA & NMDS | 2,172.502 | 568.188 | 1 | 0.583 | – | 1.04 | 0.181 (p = .000) | −0.019 (p = .750) | −0.141 (p = .012) | 0.004 (p = .962) |
| Angiosperms | bySpxSite | phyloPCOA | 2,255.31 | 650.996 | 1 | 0.653 | – | 0.816 | – | – | – | – |
| Angiosperms | bySpxSite | NMDS & site | 2,279.589 | 675.275 | 1 | 0.578 | 0.361 | – | 0.089 (p = .147) | −0.014 (p = .839) | −0.141 (p = .026) | −0.002 (p = .968) |
| Angiosperms | bySpxSite | site | 2,290.683 | 686.369 | 1 | 0.581 | 0.389 | – | – | – | – | – |
| Angiosperms | bySpxSite | NMDS | 2,504.622 | 900.308 | 1 | 0.958 | – | – | 0.089 (p = .085) | 0.040 (p = .471) | −0.117 (p = .022) | 0.002 (p = .960) |
| Angiosperms | bySpxSite | no Predictors | 2,540.015 | 935.701 | 1 | 1.002 | – | – | – | – | – | – |
Note
- All predictors used in the models are indicated by presence of regression coefficient; parameters not included in each model are indicated by blanks in the table. Models were run using two sets of data: (a) all taxa in the dataset, and (b) just the angiosperms in the dataset. Analyses were conducted by species and by site × species. Results include Deviance Information Criterion (DIC), differences between the top-ranked model and others (delta DIC), and the cumulative DIC weight (Sum DICw). p-values < 0.05 are highlighed in bold.
Average phylogenetic MNTD for the sites where more than one taxon occurred (N = 179) was 210.4 Ma, and average phylogenetic MPD was 232.5 Ma. Eight sites showed significant MPD and 10 showed significant MNTD (two-tailed test, p < .01, Figure S4). Forty-one sites sampled had at least one eudicot and one species from a different major clade.
4 DISCUSSION
The analyses presented here demonstrate that phylogenetic history—the evolutionary history that led to the plants we observe across the globe today—is a more powerful predictor of global variation in rates of in-stream leaf litter decomposition than site or the range of climate variables considered in this study. While this is perhaps not surprising at the species level, where we might expect intrinsic properties of species to swamp climatic conditions at the sites where individual plants are collected, it is remarkable that even data analysed at the site level suggests that phylogeny explains 2.2–5.8 times the variation explained by site (Table 4). This is only possible because sites are phylogenetically diverse: of the 179 sites in our study where two or more taxa were measured, the average phylogenetic MNTD and average phylogenetic MPD corresponded to divergence times of 105 and 116 Ma respectively, demonstrating that sites captured large swaths of phylogenetic history. In fact, 41 sites sampled had at least one eudicot and one species from a different major clade. As the major evolutionary transition in decomposition rates appears to have been on the branch leading to the eudicots, this co-occurrence of eudicots with non-eudicot lineages is key to the fact that site conditions as well as climatic conditions predicted less variance than even the relatively small proportion of phylogenetic variance captured by four phylogenetic eigenvectors (Table 3). Evolutionary history explains variation in traits that are highly relevant to explaining global litter decomposition rates.
Phylogenetic effects observed in this study may be related to variation in whole suites of plant traits that influence microbial decomposers and invertebrate litter consumers and that act across ecosystems, from the terrestrial environments where the plants evolved to the aquatic ecosystems where their leaves decompose. Leaf litter decomposability is influenced by a plant's position on the leaf economics spectrum (Díaz et al., 2016; Wright et al., 2004). In terrestrial systems, plant traits have been found to explain the decomposition rate of leaves (Cornwell et al., 2008) and wood (Hu et al., 2018) more than climatic factors. We were unable to include specific leaf traits as explanatory variables in our study due to inconsistent reporting of leaf traits in the published studies. However, plant traits are phylogenetically heritable (Cavender-Bares et al., 2018; Flores et al., 2014; Hao, Kuang, & Kang, 2015; Pearse & Hipp, 2012; Schmerler et al., 2012; Zanne et al., 2014), and plant phylogeny is thus expected to integrate across the entire suite of litter quality traits. Hence, our study implies that the traits that predict leaf decomposability track phylogenetic history more closely than climatic transitions among species. Our results demonstrate that phylogeny explains an important fraction of variance in plant traits similar to other studies where they were either correlated with a subset of climatic predictors (e.g. Willis, Ruhfel, Primack, Miller-Rushing, & Davis, 2008) or incompletely explained by the full set of climatic predictors expected to explain them (Li, Ives, & Waller, 2017; Pearse & Hipp, 2009).
The phylogenetic effect we find in this study is among the strongest reported to date for decomposition rates. Phylogenetic heritability of leaf decomposition in streams analysed in this study (Pagel's λ = 0.977 for 239 species globally) is higher than λ values found for the decomposability of angiosperm species studied in terrestrial environments in the UK (45 species; λ = 0.32) and central Argentina (24 species; λ = 0.70) (Díaz et al., 2013), the eastern US (λ = 0.39 for 78 leaf species, Jo, Fridley, & Frank, 2016), and Ecuador (17 species; λ = 0.80, for undamaged leaves; Cárdenas, Hättenschwiler, Valencia, Argot, & Dangles, 2015). These smaller studies had lower power to detect divergence from Brownian motion processes of trait evolution (Boettiger, Coop, & Ralph, 2012) and are biased towards reporting high λ values by the long average length of their terminal branches. Therefore, our results are particularly remarkable in the light of the fact that the high species density of our phylogenetic tree and resultant relative shortening of terminal branches increases statistical power to detect divergence from Brownian motion (compare, for example, Figure 4 of the present study with Figure 4 of Díaz et al., 2013). In fact, the phylogenetic signal we observed for in-stream decomposition is greater than the phylogenetic signals for all leaf traits (λ = 0.15–0.70 for chlorophyll, P, Ca, tannins, phenols, cellulose, N and lignin, in increasing order of λ) in a tropical forest for a recent analysis of 184 species (McManus et al., 2016) and a number of leaf traits (λ = 0.27–0.76 for leaf area, chlorophyll and leaf thickness) across 229 species in tropical and subtropical China (Yang et al., 2014).
Our large-scale synthesis demonstrates that phylogenetic variation shapes a fundamental ecosystem process across terrestrial–aquatic ecosystem boundaries, despite the fact that the species being studied evolved in one ecosystem while the process of interest—in-stream decomposition rate—is measured in another. This phylogenetic effect is detectable even at a significant distance from the source of the leaf litter and despite a broad set of environmental factors that also influence decomposition rates—water velocity, water temperature, microbial colonization and macroinvertebrate shredders. The patterns we observe here are thus not a consequence of phylogenetic niche conservatism; rather, the entire suite of phylogenetically heritable plant traits (e.g. leaf chemistry, toughness) that shape in-stream litter decomposition rates.
Our study further illustrates that the tree of life is not simply a product of adaptation to changing environments, but that evolutionary history of species and, presumably, their traits drive ecosystem processes (Cornelissen & Cornwell, 2014). The analyses presented here highlight the central role of phylogenetic history's influence on litter decomposition at a global scale and this work represents an important merger of ecosystem science with evolutionary history (Mouquet et al., 2012; Narwani, Matthews, Fox, & Venail, 2015). The result that phylogenetic signals can be observed across terrestrial–aquatic boundaries is novel and strengthens evidence for the organizing power of genes on communities and ecosystems (Jackrel & Wootton, 2014; LeRoy & Fischer, 2019; LeRoy, Whitham, Keim, & Marks, 2006; LeRoy et al., 2007; LeRoy, Wooley, & Lindroth, 2012; Whitham et al., 2006).
This study also provides further context for our continued exploration of how climate change and global warming are predicted to increase rates of aquatic decomposition through increased stream temperatures (Boyero et al., 2011; Follstad Shah et al., 2017; Tiegs et al., 2019). There are both direct and indirect influences of temperature on this key ecosystem process. Temperature can directly increase rates of decomposition through accelerated microbial growth and litter consumption; however, elevated CO2 concentrations driving temperature increases could also lower leaf litter quality via higher C:nutrient ratios and thus suppress decomposition rates (Boyero et al., 2017; Tuchman, Wetzel, Rier, Wahtera, & Teeri, 2002). Temperature can also indirectly affect decomposition through plant species’ plastic responses to climate change in terms of expressed phenotypes, adaptation in situ and shifts in species ranges (Christmas, Breed, & Lowe, 2016; Hooper et al., 2012; Kominoski et al., 2013). How these indirect changes may influence stream ecosystems is not well understood. Finally, there may be feedbacks between climate change and in-stream decomposition through carbon evasion from streams. Riverine ecosystems are important contributors to the breakdown of terrestrially derived carbon (Cole et al., 2007), and as a result, streams and rivers have been identified as major sources of carbon dioxide in the atmosphere (Raymond et al., 2013).
Streams are influenced by the biomes they are embedded in (Dodds et al., 2015; Tiegs et al., 2019), but biomes and the species assemblages within them will both shift under predicted climate change models. The selective pressures on plants that might influence litter decomposition rates (e.g. predation, defences against microbial attack, desiccation resistance) are decoupled from the physical, chemical and biological processes occurring in streams. In spite of these complexities, our study highlights the importance of riparian plant community composition on stream litter decomposition. Thus, the interactive changes in phenotypic expression, phenology and species range shifts could directly and indirectly influence global litter decomposition rates. Taken as a whole, this work highlights the significant effect of phylogenetic history on key ecosystem processes.
ACKNOWLEDGEMENTS
We thank the LTER Network for funding an initial collaboration in 2011. NAG was supported by the Department of Energy's Office of Science, Biological and Environmental Research. Oak Ridge National Laboratory is managed by UT-Battelle, LLC, for the US DOE under contract DE-AC05-00OR22725. We thank Chuck Hawkins, Sherri Johnson, Chris Swan, Amy Rosemond and Lydia Zeglin for insightful conversations about this research and for participating in the initial data compilation. We thank Bruce McCune for advice on statistical analyses. We thank Editor Dr. David Wardle and two anonymous reviewers for their comprehensive reviews of this work. All authors have confirmed that they have no relationships, financial or otherwise, that might be perceived as influencing their objectivity or that could be considered a potential source of conflict of interest.
AUTHORS’ CONTRIBUTIONS
C.J.L. and A.L.H. contributed equally to writing this paper and are joint first authors. C.J.L. conceived the study idea and coordinated collaborators. A.L.H. and K.L. compiled the phylogenetic data, and A.L.H. performed the phylogenetic comparative analyses and statistical analyses. J.J.F.S. led the initial data compilation on litter decomposition in streams, which also involved J.S.K., M.A., W.K.D., M.O.G., N.A.G., A.L., C.J.L., D.W.P.M., and J.R.W. C.J.L., A.L.H. and K.L. wrote the initial draft manuscript and all authors (C.J.L., A.L.H., K.L., J.J.F.S., J.S.K., M.A., W.K.D., M.O.G., N.A.G., A.L., D.W.P.M., R.L.S. and J.R.W.) contributed critically to subsequent drafts and gave final approval for publication. C.J.L. coordinated inputs and edits by all co-authors and reviewers.
Open Research
DATA AVAILABILITY STATEMENT
All raw and processed data, R-scripts, and figures are available at: https://github.com/andrew-hipp/decomposition-phylogeny-2019.
