Biodiversity mediates the effects of stressors but not nutrients on litter decomposition.

Understanding the consequences of ongoing biodiversity changes for ecosystems is a pressing challenge. Controlled biodiversity-ecosystem function experiments with random biodiversity loss scenarios have demonstrated that more diverse communities usually provide higher levels of ecosystem functioning. However, it is not clear if these results predict the ecosystem consequences of environmental changes that cause non-random alterations in biodiversity and community composition. We synthesized 69 independent studies reporting 660 observations of the impacts of two pervasive drivers of global change (chemical stressors and nutrient enrichment) on animal and microbial decomposer diversity and litter decomposition. Using meta-analysis and structural equation modeling, we show that declines in decomposer diversity and abundance explain reduced litter decomposition in response to stressors but not to nutrients. While chemical stressors generally reduced biodiversity and ecosystem functioning, detrimental effects of nutrients occurred only at high levels of nutrient inputs. Thus, more intense environmental change does not always result in stronger responses, illustrating the complexity of ecosystem consequences of biodiversity change. Overall, these findings provide strong evidence that the consequences of observed biodiversity change for ecosystems depend on the kind of environmental change, and are especially significant when human activities decrease biodiversity.

Introduction

n class="Species">Human activities cause global enpan>vironmenpan>tal chanpan>ges with importanpan>t consequenpan>ces for biodiversity anpan>d the funpan>ctioninpan>g of ecosystems. Understanpan>dinpan>g these consequenpan>ces is crucial for better policy anpan>d conservation strategies, which will ultimately promote pan> class="Species">human well-being too (IPBES, 2019). A key question is to what extent changes in ecosystem functioning are mediated by changes at which dimensions of biodiversity. Extensive research has demonstrated that biodiversity is needed for the stable provenance and enhancement of ecosystem processes and functions (Cardinale et al., 2012; Schuldt et al., 2018; Tilman et al., 2012). However, this body of evidence is mostly based on experiments comparing ecosystem functioning in artificial communities with varying number of species. Such experiments might not capture the complex ways by which shifts in biodiversity induced by global change ultimately affect ecosystem functioning (De Laender et al., 2016; Eisenhauer et al., 2019b).

Early biodiversity-ecosystem function (BEF) experiments typically controlled for environmental gradients, thus not accounting for the underlying drivers of biodiversity change (De Laender et al., 2016; Srivastava and Vellend, 2005; Wardle, 2016). These early experiments also focused on species richness as the sole biodiversity index, and manipulated it directly and randomly. However, environmental change will often elicit non-random changes in several facets of biodiversity (Eisenhauer et al., 2016; Giling et al., 2019; van der Plas, 2019) (community composition and population densities (Glassman et al., 2018; Spaak et al., 2017), functional diversity (Cadotte et al., 2011; Craven et al., 2018; Heemsbergen et al., 2004), trophic diversity (Soliveres et al., 2016; Wang and Brose, 2018; Zhao et al., 2019). The selective effects of environmental change en class="Gene">merge because organpan>isms differ inpan> their response to enpan>vironmenpan>tal chanpan>ge. For example, larger organpan>isms anpan>d predators are oftenpan> more negatively affected thanpan> smaller organpan>isms at lower trophic levels (Hinpan>es et al., 2015; Sheridanpan> anpan>d Bickford, 2011; Srivastava anpan>d Vellenpan>d, 2005; Voigt et al., 2007). Usinpan>g realistic extinpan>ction scenpan>arios, experimenpan>ts founpan>d contrastinpan>g effects of non-ranpan>dom shifts inpan> biodiversity on ecosystem funpan>ctioninpan>g (e.g. Cárdenpan>as et al., 2017; Jonsson et al., 2002; Melguizo-Ruiz et al., 2020; Oliveira et al., 2019; Smith anpan>d Knpan>app, 2003, Zavaleta anpan>d Hulvey, 2004). In addition, several variables that are not directly related to biodiversity control ecosystem funpan>ctions (e.g. physiological rates [Dib et al., 2020; Thakur et al., 2018] anpan>d alterations of physical anpan>d chemical conditions [De Laenpan>der et al., 2016; Gilinpan>g et al., 2019]). Whenpan> enpan>vironmenpan>tal chanpan>ge affects these mechanpan>isms, teasinpan>g out the relative importanpan>ce of biodiversity-mediated effects is complicated evenpan> more. Givenpan> the number of differenpan>t potenpan>tial mechanpan>isms, quanpan>tifyinpan>g the extenpan>t to which shifts inpan> biodiversity unpan>derpinpan> the effect of enpan>vironmenpan>tal chanpan>ge on ecosystem funpan>ctioninpan>g unpan>der real-world scenpan>arios of global chanpan>ge is a key challenpan>ge for ecology (De Laenpan>der et al., 2016; Duffy et al., 2017; Eisenpan>hauer et al., 2019b; Srivastava anpan>d Vellenpan>d, 2005; vanpan> der Plas, 2019; Wardle, 2016). Incorporatinpan>g the impacts of enpan>vironmenpan>tal chanpan>ge drivers inpan>to BEF studies anpan>d meta-anpan>alyses is anpan> importanpan>t step forward to address such questions (De Laenpan>der et al., 2016; Eisenpan>hauer et al., 2019b).

The vast majority of BEF experiments has focused on plant richness and ecosystem functions such as biomass production (van der Plas, 2019). However, litter decomposition has a tremendous importance in ecosystems and biogeochemical cycles (Follstad Shah et al., 2017). Small changes in the rate of this process can have important consequences for the overall n class="Chemical">carbon balanpan>ce. Indeed, inpan>creases inpan> decomposition rates could have positive feedback effects on climate warminpan>g by enpan>hanpan>cinpan>g C losses (Kirschbaum, 2000). The diversity of decomposers (inpan>vertebrates anpan>d microorganpan>isms that fragmenpan>t anpan>d decompose organpan>ic matter inpan> both aquatic anpan>d terrestrial systems) is crucial for litter decomposition (Eisenpan>hauer et al., 2012; García-Palacios et al., 2013; Gessnpan>er et al., 2010; Hanpan>da et al., 2014; Hättenpan>schwiler et al., 2005) anpan>d for other ecosystem funpan>ctions as well (Eisenpan>hauer et al., 2019a; Lefcheck et al., 2015; Schuldt et al., 2018). Despite the importanpan>ce of decomposers, BEF experimenpan>ts focusinpan>g on litter decomposition more oftenpan> addressed the inpan>fluenpan>ce of planpan>t litter diversity thanpan> of decomposers (Gessnpan>er et al., 2010; Toninpan> et al., 2018). In a meta-anpan>alysis, decomposer diversity had a greater effect on decomposition thanpan> the diversity of planpan>t litter (Srivastava et al., 2009), although also weak anpan>d neutral effects have beenpan> reported (vanpan> der Plas, 2019). Facilitation anpan>d complemenpan>tarity through niche partitioninpan>g are primary mechanpan>isms unpan>derlyinpan>g the positive relationship betweenpan> decomposer diversity anpan>d decomposition (Gessnpan>er et al., 2010; Hättenpan>schwiler et al., 2005; Toninpan> et al., 2018). Experimenpan>ts conducted inpan> natural conditions anpan>d reflectinpan>g realistic extinpan>ction scenpan>arios are still relatively scarce, anpan>d demonstrate contrastinpan>g effects of non-ranpan>dom shifts inpan> decomposer diversity on decomposition (Cárdenpan>as et al., 2017; Jonsson et al., 2002; Melguizo-Ruiz et al., 2020; Wenpan>isch et al., 2017). The need to quanpan>tify enpan>vironmenpan>tal chanpan>ge effects on decomposer diversity, along with potenpan>tial knpan>ock-on effects on litter decomposition, is therefore particularly pressinpan>g.

There is a variety of environmental change drivers, and different types of drivers may have diverse effects on biodiversity and ecosystem functions (De Laender et al., 2016; Dib et al., 2020). We postulate that there are two main categories of environmental change: stressors anpan>d resource shifts. While pan> class="Disease">stressors cannot be consumed, and act as conditions that alter growth rates (e.g., temperature, drought, chemical stressors), resources are by definition consumed (e.g., CO2 or mineral nutrients), which has important implications for how they should enter theory (Chase and Leibold, 2003; De Laender, 2018). Chemical stressors and nutrient enrichment are important case studies of environmental stressors and resource enrichment, because their presence is increasing rapidly (Bernhardt et al., 2017) and they are projected to have severe effects on biodiversity (Mazor et al., 2018). They are also of particular relevance for decomposer communities. Chemical stressors such as metals and pesticides decrease the diversity, abundance, growth and activity of decomposers across terrestrial and aquatic systems (e.g. Hogsden and Harding, 2012; Pelosi et al., 2014; Schäfer, 2019). In contrast, nutrient enrichment can have positive impacts on the abundance and physiological rates of decomposer organisms by reducing resource limitations (Treseder, 2008), but at the same time decrease decomposer diversity (Lecerf and Chauvet, 2008; Woodward et al., 2012). Across ecosystems, stressors and nutrients can exert opposite impacts on litter decomposition rates, with decreases in response to chemical stressors but increases following nutrient enrichment (Ferreira et al., 2015; Ferreira et al., 2016). In addition, decomposition involves both microorganisms and invertebrates (Bardgett and van der Putten, 2014; Gessner et al., 2010; Hättenschwiler et al., 2005) that may respond differently to stressors and nutrients with a higher sensitivity of invertebrates than microorganisms (Peters et al., 2013; Siebert et al., 2019). Although many published case studies report shifts in decomposer diversity and in rates of litter decomposition at sites impacted by stressors and nutrients, biodiversity-mediated effects have not yet been quantified across systems.

Here we addressed the question if the effects of stressors anpan>d pan> class="Gene">nutrient enrichment on decomposer diversity and abundance explain the response of litter decomposition to these two types of pervasive environmental change drivers (Figure 1). We synthesized 69 published case studies reporting the impact of stressors (metals, pesticides) and nutrients (nitrogen or phosphorous additions) on litter decomposition and on decomposer diversity (taxa richness, Shannon diversity, evenness) or abundance (density, biomass) at sites differing in stressor or nutrient levels. Our comprehensive global dataset of 660 observations encompasses studies across taxonomic groups (animal (soil micro-, meso- and macrofauna, stream macroinvertebrates) and microbial (fungi and bacteria) decomposers), ecosystems (aquatic and terrestrial), and study types (experimental and observational) (Figure 2). We quantified the effect size of environmental change on decomposer diversity or abundance and on litter decomposition within studies using correlation coefficients between stressor or nutrient levels and decomposer diversity, abundance, and litter decomposition. We also characterized stressor and nutrient intensities, and standardized their levels in water, soil, or sediment using environmental quality criteria issued by environmental authorities (e.g. ECHA, USEPA, UKTAG). Using meta-analysis and structural equation modelling (SEM), we first compared the overall effects of stressors and nutrients on decomposers and decomposition across systems and studies (first meta-analysis), and second, addressed to what extent changes in decomposer diversity and abundance mediate the impacts of these two contrasting drivers of environmental change on decomposition (second meta-analysis and SEM). Third, we explored the effects of three main moderators on decomposers diversity, abundance, and decomposition responses, as found in the second meta-analysis: stressor or nutrient intensity, taxonomic group (animal vs. microbes) and study type (experimental vs. observational studies).

Figure 1.

Schematic representation of the structural hypotheses tested in this study.

Green arrows depict expected positive effects, red arrows represent negative effects. Stressors and nutrients are hypothesized to decrease decomposer diversity. The response of decomposers diversity to environmental change drivers determines the response of decomposition (Srivastava et al., 2009). Nutrients are hypothesized to increase decomposer abundance. Stressors and nutrients can affect litter decomposition independent of changes in decomposer diversity and abundance, especially through changes in physiological activity (De Laender et al., 2016; Giling et al., 2019).

Figure 2.

Description of the data used in the present meta-analysis.

(A) Countries represented and corresponding number of observations, (B) decomposer diversity and abundance metrics covered, and (C) ecosystem types and decomposer taxonomic groups (animals: soil micro-, meso-, macro-fauna, stream macroinvertebrates; and microbial decomposers: fungi and bacteria) represented.

Green arrows depict expected positive effects, red arrows represent negative effects. n class="Disease">Stressors anpan>d pan> class="Gene">nutrients are hypothesized to decrease decomposer diversity. The response of decomposers diversity to environmental change drivers determines the response of decomposition (Srivastava et al., 2009). Nutrients are hypothesized to increase decomposer abundance. Stressors and nutrients can affect litter decomposition independent of changes in decomposer diversity and abundance, especially through changes in physiological activity (De Laender et al., 2016; Giling et al., 2019).

We expected that chemical stressors anpan>d pan> class="Gene">nutrients would have contrasting effects on decomposer diversity and abundance, and on litter decomposition across ecosystems and studies (Figure 1). We hypothesized that chemical stressors generally decrease decomposer diversity, abundance (Hogsden and Harding, 2012; Petrin et al., 2008), and litter decomposition rates (Ferreira et al., 2016; Peters et al., 2013), and that nutrients generally decrease decomposer diversity (Lecerf and Chauvet, 2008; Woodward et al., 2012) but increase decomposer abundance and litter decomposition rates (based on physiological effects and decreasing resource limitations (Bergfur et al., 2007; Ferreira et al., 2015; Treseder, 2008; Woodward et al., 2012). We further hypothesized that litter decomposition responses to environmental change depend on changes in decomposer diversity and abundance, and expected an overall positive relationship independent of environmental change intensity (Srivastava et al., 2009).

Results

Description of the data and overall patterns

The final dataset contained 69 (case) studies from 59 publications, representing 660 observations. Data were mostly from Europe (44 ; 443 (studies; observations)) and n class="Chemical">North anpan>d South Apan> class="Gene">merica (19; 168), while Asia (2; 9) and Oceania (4; 40) were less well represented (Figure 2A). The studies covered aquatic (55; 388) and terrestrial systems (14; 272) (Figure 2C), and used observational (43; 336) or experimental approaches (26; 324). Studies reported abundance (66; 463) or diversity responses (48; 197) (Figure 2B) of soil and benthic invertebrates (48; 509) and microbes (fungi and bacteria) associated with litter materials (36; 151) (Figure 2C). Chemical stressors were mostly metals (13; 257) and pesticides (12; 66) associated with industrial activities, accidental spills, and agricultural practices. Nutrient enrichment studies addressed fertilization by various N and/or P forms (26; 175), and eutrophication due to agricultural runoffs (10; 59) or wastewater effluents (4; 44). There was no study reporting nutrient enrichment impacts on soil decomposer diversity in the dataset. Funnel plots and intercepts of Egger’s regression showed evidence for positive publication bias in nutrient enrichment studies reporting decomposer abundance (Appendix 2—figure 1; Appendix 2—figure 2; Appendix 2—table 1). No publication bias was detected in the other datasets.

Appendix 2—figure 1.

Assessment of publication bias.

Stressors: Funnel plots of each response variables (decomposer diversity, abundance and decomposition) in the two datasets (stressors - diversity and stressors - abundance). Meta-analytic models included the effect of stressor intensity (standardized levels) as a covariate.

Appendix 2—figure 2.

Assessment of publication bias.

Nutrients: Funnel plots of each response variables (decomposer diversity, abundance and decomposition) in the two datasets (stressors - diversity and stressors - abundance). Meta-analytic models included the effect of nutrient intensity (standardized levels) as a covariate.

Appendix 2—table 1.

Assessment of publication bias.

Results from Egger’s regressions showing the intercepts, standard error (SE) and p-value of regressions between standard normal deviate of each response variable (effect sizes) and the inverse of their standard errors. Models also included stressor or nutrient intensity as a covariate.

Dataset	Variable	Publication bias p	Publication bias	Intercept	SE
Stressors - Biodiv	Biodiversity	0.10	no	−1.36	0.83
Stressors - Biodiv	Decomposition	0.58	no	−1.07	1.94
Stressors - Abdc	Abundance	0.14	no	−1.49	1.02
Stressors - Abdc	Decomposition	0.68	no	−0.67	1.60
Nutrients - Biodiv	Biodiversity	0.37	no	0.76	0.86
Nutrients - Biodiv	Decomposition	0.19	no	3.35	2.55
Nutrients - Abdc	Abundance	0.08	no	1.21	0.70
Nutrients - Abdc	Decomposition	<0.001	pub. bias	5.31	1.45

We found largely contrasting effects of stressors anpan>d pan> class="Gene">nutrients on each of the three response variables in a first-level meta-analysis comparing the overall effects of the two drivers of environmental change (Figure 3, Appendix 2—table 2). Chemical stressors overall decreased decomposer diversity, abundance and litter decomposition across studies (Figure 3). Nutrient enrichment tended to decrease decomposer diversity but to increase abundance, and decomposition, although these trends were not significant as indicated by confidence intervals of the grand mean effects overlapping with zero (Figure 3).

Figure 3.

Grand mean effect sizes of chemical stressors and nutrient enrichment on decomposer diversity (taxa richness and diversity indices), abundance (density and biomass), and litter decomposition.

Effect sizes are z-transformed correlation coefficients. Error bars show 95% confidence intervals. Numbers in parentheses indicate number of studies and observations, respectively. Symbols show the significance level for the comparison between mean effect size and zero (***p<0.001; *p<0.05). For full model results, see Appendix 2—table 2.

Appendix 2—table 2.

First level meta-analysis comparing the effects of chemical stressors and nutrient enrichment.

Results of Wald-type chi-square tests comparing the grand mean effect sizes of the three response variables (decomposer diversity, abundance and litter decomposition) between chemical stressors and nutrient enrichment.

Response	QM	Df	N	p-value
Diversity	25.65	2	174	<0.001
Abundance	7.92	2	424	0.019
Litter decomposition	17.61	2	165	<0.001

Effect sizes are z-transformed correlation coefficients. Error bars show 95% confidence intervals. n class="Chemical">Numbers inpan> parentheses inpan>dicate number of studies anpan>d observations, respectively. Symbols show the significanpan>ce level for the comparison between meanpan> effect size anpan>d zero (***p<0.001; *p<0.05). For full model results, see Appendix 2—table 2.

Biodiversity-mediated effects of stressors and nutrients on litter decomposition

The responses of decomposition and of decomposer diversity and abundance to chemical n class="Disease">stressors were correlated: decreases inpan> decomposition were associated with decreases inpan> decomposer diversity anpan>d abunpan>danpan>ce (Figure 4 upper panpan>els). We did not finpan>d such a relationship for pan> class="Gene">nutrients. Instead, a range of positive and negative responses of decomposer diversity, abundance, and decomposition to nutrients were found, without significant associations between them (Figure 4 lower panels). In addition, when decomposer diversity and abundance responses to nutrients were close to zero, there was a wide range of decomposition responses (intercepts from Figure 4 lower panels).

Figure 4.

Relationship between the responses of decomposition and decomposer diversity and abundance to chemical stressors and nutrient enrichment.

Variables are effect sizes (z-transformed correlation coefficients) of stressors or nutrients on litter decomposition and on animal and microbial decomposer diversity (left panels) or abundance/biomass (right panels). Gray symbols are individual observations of effect sizes; Colored symbols indicate mean effect size on diversity or abundance across individual observations for a unique litter decomposition measurement (see methods). Lines represent meta-regressions between effect sizes for decomposition and decomposers, where solid lines are statistically significant (p<0.05), dashed lines are non-significant (p>0.05), and thin lines depict the regression's confidence interval. QM and p represent the model heterogeneity and p-values of the meta-regressions, respectively, with sample size (number of studies; number of observations).

Variables are effect sizes (z-transformed correlation coefficients) of n class="Disease">stressors or pan> class="Gene">nutrients on litter decomposition and on animal and microbial decomposer diversity (left panels) or abundance/biomass (right panels). Gray symbols are individual observations of effect sizes; Colored symbols indicate mean effect size on diversity or abundance across individual observations for a unique litter decomposition measurement (see methods). Lines represent meta-regressions between effect sizes for decomposition and decomposers, where solid lines are statistically significant (p<0.05), dashed lines are non-significant (p>0.05), and thin lines depict the regression's confidence interval. QM and p represent the model heterogeneity and p-values of the meta-regressions, respectively, with sample size (number of studies; number of observations).

According to our overarching hypothesis, the SEM indicated that the effects of stressors on litter decomposition were mediated by shifts inpan> decomposer diversity anpan>d abunpan>danpan>ce. Includinpan>g the direct paths from decomposer diversity or abunpan>danpan>ce to litter decomposition improved both the models accordinpan>g to mediation tests anpan>d AIC comparisons (Figure 5). In addition, the path coefficienpan>ts from diversity anpan>d abunpan>danpan>ce to the decomposition response to pan> class="Disease">stressors had (standardized) values higher than 0.1 (Figure 5) and were statistically different from zero (Appendix 2—table 3). However, in contrast to chemical stressors, the SEM did not support biodiversity-mediated effects of nutrient enrichment on litter decomposition. While the mediation test and AIC indicated that the decomposer diversity-mediated path improved the model (Figure 5), the path coefficient was not significantly different from 0 (Appendix 2—table 3). The decomposer abundance-mediated path of nutrients was not supported by the data: an SEM without the direct path from decomposer abundance to decomposition could not be rejected based on the mediation test (Figure 5), and including this path did not improve the model according to the AIC comparison. Besides, we found publication bias in this dataset (Appendix 2—figure 2, Appendix 2—table 1), and model check indicated that the residuals of the nutrients-abundance model were non-independent from the fitted values. Thus, the results from this model are reported here for comparison purposes only.

Figure 5.

Decomposer diversity and abundance explained litter decomposition response to chemical stressors but not to nutrient enrichment.

Structural equation models investigating decomposer diversity- or abundance-mediated effects of chemical stressors and nutrient enrichment on litter decomposition across 69 studies. Arrows represent relationships between stressor or nutrient intensity levels, and effect sizes of stressors or nutrients on litter decomposition and on decomposer diversity (taxa richness, Shannon diversity, or evenness: left panels) or abundance and biomass (right panels). Values along the arrows are standardized path coefficients. Green, red, and gray arrows indicate positive, negative, and non-significant relationships, respectively. Curved arrows depict the indirect effects of stressors or nutrients on decomposition as mediated by diversity or abundance. Arrow widths are scaled relative to the magnitude of standardized path coefficients. C statistic, P-value (P<0.05 indicate poor model fit), and sample sizes (number of studies; number of observations). Results of mediation tests: comparison with models omitting the path from diversity or abundance to decomposition (ΔAIC < −2 indicates that reduced models were not consistent with the data).

Appendix 2—table 3.

Summary table of structural equation modelling (SEM) analysis.

Unstandardized path coefficients from SEMs for the four datasets: Stressors - Biodiversity (Biodiv), Stressors - Abundance (Abdc), Nutrients - Biodiversity and Nutrients, Abundance. SEMs also incorporated categorical predictors (study type, taxonomic group and diversity metric, see Materials and methods).

Dataset	Response	Predictor	Estimate	SE	Crit.value	Df	p-Value
Stressors - Biodiv	Decomposition	Diversity	0.42	0.17	2.50	19	0.022
Stressors - Biodiv	Decomposition	Stressor intensity	−0.02	0.04	−0.47	19	0.643
Stressors - Biodiv	Diversity	Stressor intensity	−0.10	0.04	−2.44	18	0.025
Stressors - Abdc	Decomposition	Abundance	0.24	0.08	2.97	25	0.007
Stressors - Abdc	Decomposition	Stressor intensity	−0.01	0.03	−0.41	25	0.683
Stressors - Abdc	Abundance	Stressor intensity	0.00	0.05	0.03	25	0.977
Nutrients - Biodiv	Decomposition	Diversity	0.01	0.11	0.06	20	0.951
Nutrients - Biodiv	Decomposition	Nutrient intensity	−0.08	0.06	−1.21	20	0.239
Nutrients - Biodiv	Diversity	Nutrient intensity	−0.25	0.07	−3.51	19	0.002
Nutrients - Abdc	Decomposition	Abundance	0.08	0.10	0.76	44	0.451
Nutrients - Abdc	Decomposition	Nutrient intensity	−0.12	0.05	−2.16	44	0.037
Nutrients - Abdc	Abundance	Nutrient intensity	−0.06	0.06	−1.00	44	0.321

Structural equation models investigating decomposer diversity- or abundance-mediated effects of chemical stressors anpan>d pan> class="Gene">nutrient enrichment on litter decomposition across 69 studies. Arrows represent relationships between stressor or nutrient intensity levels, and effect sizes of stressors or nutrients on litter decomposition and on decomposer diversity (taxa richness, Shannon diversity, or evenness: left panels) or abundance and biomass (right panels). Values along the arrows are standardized path coefficients. Green, red, and gray arrows indicate positive, negative, and non-significant relationships, respectively. Curved arrows depict the indirect effects of stressors or nutrients on decomposition as mediated by diversity or abundance. Arrow widths are scaled relative to the magnitude of standardized path coefficients. C statistic, P-value (P<0.05 indicate poor model fit), and sample sizes (number of studies; number of observations). Results of mediation tests: comparison with models omitting the path from diversity or abundance to decomposition (ΔAIC < −2 indicates that reduced models were not consistent with the data).

The magnitude of the biodiversity-mediated effects of chemical stressors on decomposition was stronger thanpan> that of the direct effects of pan> class="Disease">stressor intensity on decomposition. The indirect effect of stressors on decomposition mediated by diversity (i.e. mathematical product of the standardized paths from stressor intensity to decomposer diversity and from diversity to decomposition Figure 5) was higher than the direct effect of stressors on decomposition, while the abundance-mediated effect of stressors was negligible (Figure 5). In the case of nutrient enrichment, however, decomposition responses were not explained by shifts in decomposer diversity and abundance, and the direct effects of nutrient intensity dominated the total effect (Figure 5). Finally, between-model comparisons (based on unstandardized path coefficients [Grace, 2006]) revealed that decomposer diversity was a stronger driver of decomposition response to stressors than decomposer abundance (unstandardized paths were 0.42 and 0.24 for diversity and abundance, respectively, Appendix 2—table 3).

Sensitivity analyses revealed that the results were robust to the inclusion of approximated standard deviations (Appendix 3—table 1; Appendix 3—table 2), and extreme values of effect sizes (Appendix 3—table 3; Appendix 3—table 4). We found partially different results when using log-response ratios as effect sizes (Appendix 3—table 5; Appendix 3—table 6), due to lower sample sizes and en class="Gene">mergenpan>ce of extreme values inpan> these datasets. In addition, the log-response ratio is probably senpan>sitive to the various metrics of biodiversity, abunpan>danpan>ce, anpan>d decomposition covered by the inpan>dividual studies that we inpan>cluded, while correlation coefficienpan>ts better accommodate such discrepanpan>cies (Koricheva et al., 2013).

Appendix 3—table 1.

Results of mediation tests from structural equation modeling (SEM) analysis based on data without approximated standard deviations.

C statistic and associated p-value for SEM without the path from biodiversity or abundance to decomposition for the four datasets: Stressors - Diversity, Stressors - Abundance, Nutrients - Diversity and Nutrients - Abundance. ΔAIC is the difference in AIC score between models with and without biodiversity- or abundance-mediated effects.

Dataset	C statistic	Df	p-value	ΔAIC	No. of studies	N
Stressors, Biodiv	12.42	6	0.053	−8.32	16	58
Stressors, Abdc	10.15	4	0.038	−6.82	23	216
Nutrient, Biodiv	13.33	6	0.038	−1.46	21	67
Nutrient, Abdc	3.82	4	0.432	−0.12	32	127

Appendix 3—table 2.

Summary table of structural equation modeling (SEM) analysis based on data without approximated standard deviations.

Standardized (Std.est.) and unstandardized estimate (Est.) path coefficients from SEMs for the four datasets.

Dataset	Response	Predictor	Std.est.	Est.	SE	Crit.value	Df	p-value
Stress., Biodiv	Decomposition	Diversity	0.52	0.50	0.16	3.16	12	0.008
Stres., Biodiv	Decomposition	Stressor intensity	−0.26	−0.05	0.03	−1.54	12	0.148
Stress., Biodiv	Diversity	Stressor intensity	−0.39	−0.08	0.04	−1.89	11	0.085
Stress., Abdc	Decomposition	Abundance	0.40	0.27	0.09	2.91	19	0.009
Stress., Abdc	Decomposition	Stressor intensity	−0.11	−0.02	0.03	−0.77	19	0.450
Stress., Abdc	Abundance	Stressor intensity	0.08	0.03	0.06	0.46	19	0.649
Nut., Biodiv	Decomposition	Diversity	−0.04	−0.04	0.12	−0.35	10	0.732
Nut., Biodiv	Decomposition	Nutrient intensity	−0.31	−0.14	0.09	−1.52	10	0.161
Nut., Biodiv	Diversity	Nutrient intensity	−0.49	−0.23	0.10	−2.39	9	0.040
Nut., Abdc	Decomposition	Abundance	0.05	0.04	0.13	0.33	29	0.742
Nut., Abdc	Decomposition	Nutrient intensity	−0.26	−0.12	0.06	−1.91	29	0.066
Nut., Abdc	Abundance	Nutrient intensity	−0.20	−0.10	0.07	−1.40	29	0.173

Appendix 3—table 3.

Results of mediation tests from structural equation modeling (SEM) analysis based on data excluding extreme values of effect sizes.

C statistic and associated p-value for SEM without the path from biodiversity or abundance to decomposition for the four datasets: Stressors - Diversity, Stressors - Abundance, Nutrients - Diversity and Nutrients - Abundance. ΔAIC is the difference in AIC score between models with and without biodiversity- or abundance-mediated effects.

Dataset	C statistic	Df	p-value	ΔAIC	No. of studies	N
Stressors, Biodiv	10.18	6	0.117	−6.71	22	94
Stressors, Abdc	7.39	4	0.117	−4.23	27	254
Nutrient, Biodiv	14.80	6	0.022	−4.85	26	93
Nutrient, Abdc	2.74	4	0.603	0.15	35	159

Appendix 3—table 4.

Summary table of structural equation modelling (SEM) analysis based on data excluding extreme values of effect sizes.

Standardized (Std.est.) and unstandardized estimate (Est.) path coefficients from SEMs for the four datasets.

Dataset	Response	Predictor	Std.est.	Est.	SE	Crit.value	Df	p-value
Stress., Biodiv	Decomposition	Diversity	0.41	0.40	0.18	2.20	18	0.041
Stress., Biodiv	Decomposition	Stressor intensity	−0.04	−0.01	0.04	−0.24	18	0.814
Stress., Biodiv	Diversity	Stressor intensity	−0.44	−0.10	0.04	−2.75	17	0.014
Stress., Abdc	Decomposition	Abundance	0.30	0.24	0.11	2.24	23	0.035
Stress., Abdc	Decomposition	Stressor intensity	0.05	0.01	0.03	0.35	23	0.731
Stress., Abdc	Abundance	Stressor intensity	0.00	0.00	0.04	−0.02	23	0.980
Nut., Biodiv	Decomposition	Diversity	0.00	0.00	0.11	0.02	19	0.986
Nut., Biodiv	Decomposition	Nutrient intensity	−0.18	−0.08	0.06	−1.30	19	0.210
Nut., Biodiv	Diversity	Nutrient intensity	−0.53	−0.24	0.07	−3.36	18	0.003
Nut., Abdc	Decomposition	Abundance	0.00	0.00	0.09	0.04	37	0.968
Nut., Abdc	Decomposition	Nutrient intensity	−0.38	−0.13	0.04	−3.26	37	0.002
Nut., Abdc	Abundance	Nutrient intensity	−0.24	−0.09	0.05	−1.73	37	0.092

Appendix 3—table 5.

Results of mediation tests from structural equation modeling (SEM) analysis based on data using log-response ratio as an effect size.

C statistic and associated p-value for SEM without the path from biodiversity or abundance to decomposition for the four datasets: Stressors - Diversity, Stressors - Abundance, Nutrients - Diversity and Nutrients - Abundance. ΔAIC is the difference in AIC score between models with and without biodiversity- or abundance-mediated effects.

Dataset	C statistic	Df	p-value	ΔAIC	No. of studies	N
Stressors, Biodiv	4.11	6	0.662	−0.02	22	70
Stressors, Abdc	5.59	4	0.232	−2.22	37	150
Nutrient, Biodiv	8.03	6	0.236	−2.08	14	78
Nutrient, Abdc	3.41	4	0.492	−0.44	21	307

Appendix 3—table 6.

Summary table of structural equation modeling (SEM) analysis based on data using log-response ratio as an effect size.

Standardized (Std.est.) and unstandardized estimate (Est.) path coefficients from SEMs for the four datasets.

Dataset	Response	Predictor	Std.est	Est.	SE	Crit.value	Df	p-value
Stress., Biodiv	Decomposition	Diversity	0.18	0.12	0.15	0.80	15	0.437
Stress., Biodiv	Decomposition	Stressor intensity	−0.24	−0.05	0.04	−1.47	15	0.163
Stress., Biodiv	Diversity	Stressor intensity	−0.35	−0.12	0.03	−4.17	15	0.001
Stress., Abdc	Decomposition	Abundance	0.14	0.04	0.05	0.86	28	0.396
Stress., Abdc	Decomposition	Stressor intensity	0.09	0.02	0.04	0.55	28	0.586
Stress., Abdc	Abundance	Stressor intensity	−0.14	−0.11	0.11	−1.03	28	0.312
Nut., Biodiv	Decomposition	Diversity	0.29	0.19	0.10	1.80	14	0.094
Nut., Biodiv	Decomposition	Nutrient intensity	−0.15	−0.07	0.08	−0.96	14	0.352
Nut., Biodiv	Diversity	Nutrient intensity	−0.20	−0.16	0.07	−2.11	14	0.054
Nut., Abdc	Decomposition	Abundance	0.06	0.04	0.06	0.59	42	0.559
Nut., Abdc	Decomposition	Nutrient intensity	−0.36	−0.16	0.05	−3.08	42	0.004
Nut., Abdc	Abundance	Nutrient intensity	−0.01	0.00	0.08	−0.08	42	0.935

Response of animal and microbial decomposers and decomposition to stressor and nutrient intensity

Despite the overall negative effects of stressors on decomposition, negative responses inpan> decomposition were not associated with higher pan> class="Disease">stressor intensity (Figure 5, Figure 6). This result held for two complementary approaches: multivariate SEM (Figure 5) that relied on data resampling to account for replicated values of decomposition matching several decomposer responses (e.g. for different taxa in the same litterbag), and meta-regressions (Figure 6) where data resampling was not necessary (see Materials and methods). There was mixed support for a stressor intensity effect on decomposer diversity across the two approaches: decomposer diversity responses decreased with stressor intensity according to the SEM (Figure 5), but this trend was not significant according to the second level meta-analysis (Figure 6). Similar slopes were obtained both with the SEM relying on data resampling (the slope of the relationship was −0.10 ± 0.04, Appendix 2—table 3) and with the meta-regression (the slope was −0.05 ± 0.03). The differences between the two approaches can be explained by the different data included. Decomposer abundance responses were not associated to stressor intensity in both the SEM and meta-regression approaches (Figure 5, Figure 6). We found different patterns for nutrient enrichment, where decomposition responses decreased with nutrient intensity (Figure 5, Figure 6), from positive effects at low intensity to negative effects at higher intensity (Figure 6). A similar pattern was observed for decomposer diversity, where responses decreased with nutrient intensity from positive to neutral to negative responses at high nutrient levels (Figure 6). Nutrient intensity, however, did not explain the responses of decomposer abundance (Figure 5, Figure 6), and both positive and negative responses were found at high nutrient levels.

Figure 6.

Decomposer and decomposition responses to the intensity levels of chemical stressors and nutrient enrichment.

Values are effect sizes (z-transformed correlation coefficients). Stressor or nutrient intensity represents the standardized level of environmental change in the treatment with the highest level (values < 0: observed level below quality criteria considered to be safe for the environment; values > 0: observed level above quality criteria). Point size is proportional to the inverse of the variance in effect size. Lines are the slopes and 95% confidence intervals from bivariate meta-regressions, with associated QM statistics, p-value and sample size (number of studies; number of observations).

Values are effect sizes (z-transformed correlation coefficients). n class="Disease">Stressor or pan> class="Gene">nutrient intensity represents the standardized level of environmental change in the treatment with the highest level (values < 0: observed level below quality criteria considered to be safe for the environment; values > 0: observed level above quality criteria). Point size is proportional to the inverse of the variance in effect size. Lines are the slopes and 95% confidence intervals from bivariate meta-regressions, with associated QM statistics, p-value and sample size (number of studies; number of observations).

The meta-analysis further revealed clear discrepancies between the response of animal and microbial (fungi and bacteria) decomposers to stressors anpan>d pan> class="Gene">nutrients. Animal decomposers responded more strongly to chemical stressors than microbial decomposers. The mean effects of chemical stressors on animal decomposer diversity and abundance were more negative than that on microbial decomposers, confirmed by Wald type tests of the second-level meta-analyses (Figure 7 upper panels, Appendix 2—table 4). Animal decomposers overall decreased in diversity but increased in abundance in response to nutrient enrichment (Figure 7, lower panels). On the other hand, the mean effects of nutrients on microbial decomposer diversity and abundance had lower magnitudes compared to animals (Appendix 2—table 4), with confidence intervals overlapping with zero (Figure 7 lower left panel). Finally, there was no clear difference between observational and experimental studies (Figure 7, Appendix 2—table 4), and between biodiversity responses in terms of taxa richness or of diversity indices (Appendix 2—table 4).

Figure 7.

Moderator effects on decomposer diversity and abundance responses to chemical stressors and nutrient enrichment.

Responses of decomposer diversity (taxa richness and diversity indices) and abundance (densities and biomass) to stressors and nutrients according to the taxonomic group (animals and microbes) and study type (Expe. = experimental; Obs. = observational studies). Values are mean effect sizes (z-transformed correlation coefficients) and 95% confidence intervals derived from meta-analytic models. Sample sizes are reported for each moderator: (number of studies; number of observations).

Appendix 2—table 4.

Main effects of categorical predictors on decomposer diversity, abundance and decomposition in the four datasets: Stressors - Biodiversity (Biodiv), Stressors - Abundance (Abdc), Nutrients - Biodiversity and Nutrients, Abundance.

Results are QM statistics and associated p-values of the second-level meta-analyses.

Dataset	Response	Predictor	QM	p-value
Stressors - Biodiv	Diversity	Taxonomic group	4.80	0.028
Stressors - Abdc	Abundance	Taxonomic group	10.10	0.001
Nutrients - Biodiv	Diversity	Taxonomic group	12.77	<0.001
Nutrients - Abdc	Abundance	Taxonomic group	4.53	0.033
Stressors - Biodiv	Diversity	Study type	1.89	0.169
Stressors - Abdc	Abundance	Study type	0.92	0.338
Nutrients - Biodiv	Diversity	Study type	0.24	0.625
Nutrients - Abdc	Abundance	Study type	0.98	0.323
Stressors - Biodiv	Diversity	Diversity metric	1.67	0.196
Nutrients - Biodiv	Diversity	Diversity metric	2.35	0.125
Stressors - Biodiv	Decomposition	Study type	0.16	0.693
Stressors - Abdc	Decomposition	Study type	1.85	0.174
Nutrients - Biodiv	Decomposition	Study type	2.69	0.101
Nutrients - Abdc	Decomposition	Study type	0.18	0.674

Responses of decomposer diversity (taxa richness and diversity indices) and abundance (densities and biomass) to n class="Disease">stressors anpan>d pan> class="Gene">nutrients according to the taxonomic group (animals and microbes) and study type (Expe. = experimental; Obs. = observational studies). Values are mean effect sizes (z-transformed correlation coefficients) and 95% confidence intervals derived from meta-analytic models. Sample sizes are reported for each moderator: (number of studies; number of observations).

Discussion

The present synthesis brings new insights into how changes in decomposer biodiversity induced by two pervasive drivers of environmental change ultimately affect decomposition. We find concomitant changes in biodiversity and decomposition under the influence of chemical stressors but not pan> class="Gene">nutrient enrichment, highlighting that real-world patterns relating shifts in biodiversity and ecosystem functioning depend on the type of environmental change. In fact, we observed significant correlations between effects on biodiversity and ecosystem function in a scenario where chemical stressors caused a significant decline in biodiversity. In contrast, in cases where nutrient enrichment caused variable responses in biodiversity, relationships between biodiversity and ecosystem function responses were weaker. It remains an understudied but important question if results of controlled BEF experiments are applicable to non-random changes in biodiversity caused by human activities (e.g. De Laender et al., 2016; Duffy et al., 2017; Eisenhauer et al., 2019b; Srivastava and Vellend, 2005; van der Plas, 2019; Wardle, 2016). The present results provide strong empirical evidence for significant real-world BEF relationships when environmental changes decrease biodiversity.

Biodiversity-mediated effects of chemical stressors on decomposition

Chemical stressors caused consistenpan>t reductions inpan> decomposer diversity anpan>d abunpan>danpan>ce as well as inpan> litter decomposition rates, inpan> linpan>e with several previous case studies (Beketov et al., 2013; Malaj et al., 2014) anpan>d meta-anpan>alyses (Ferreira et al., 2016; Peters et al., 2013). Addinpan>g to the previous knpan>owledge, the presenpan>t meta-anpan>alysis shows that chanpan>ges inpan> decomposer diversity anpan>d abunpan>danpan>ce explainpan>ed the decomposition response to pan> class="Disease">stressors, providing evidence for the expectation that shifts in biodiversity mediate the impact of chemical stressors on decomposition. We acknowledge that despite the SEM analysis, the approach conducted here remains correlative. However, our study builds on a body of experimental and observational evidence that already demonstrated that more diverse and abundant decomposer communities support higher decomposition rates, albeit not under the influence of environmental change (e.g. García-Palacios et al., 2013; Handa et al., 2014).

We especially complement a previous meta-analysis showing the importance of decomposer diversity for decomposition across experiments manipulating the richness of invertebrate and microbial decomposer communities (Srivastava et al., 2009). We extend on this and show that non-random biodiversity losses induced by n class="Disease">stressors are closely associated with decreases inpan> decomposition across a wide ranpan>ge of studies. A recenpan>t review poinpan>ted out that inpan> naturally assembled terrestrial communpan>ities, studies more oftenpan> founpan>d neutral anpan>d to a lesser extenpan>t positive relationships betweenpan> decomposer diversity anpan>d decomposition (vanpan> der Plas, 2019). In that review, communpan>ities were not inpan>fluenpan>ced by enpan>vironmenpan>tal chanpan>ge drivers, anpan>d the vote counpan>tinpan>g approach used is senpan>sitive to the statistical power of inpan>dividual studies anpan>d could have inpan>creased the probability of finpan>dinpan>g non-signpan>ificanpan>t relationships (Koricheva et al., 2013). In linpan>e with our finpan>dinpan>gs, anpan> experimenpan>t mimickinpan>g the sequenpan>ce inpan> which freshpan> class="Chemical">water invertebrate decomposers are lost after disturbances showed that decreasing non-randomly the number of species decreased decomposition rates (Jonsson et al., 2002).

Biodiversity-ecosystem function experiments manipulating biodiversity directly are key to understand the mechanisms involved in this relationship (Eisenhauer et al., 2016), especially because they control for the effects of environmental heterogeneity or abundance. However, in real-world scenarios, environmental change drivers affect both biodiversity and abundance simultaneously. As demonstrated here, this is especially the case for stressors that decrease decomposer diversity anpan>d abunpan>danpan>ce (Hogsdenpan> anpan>d Hardinpan>g, 2012). The abunpan>danpan>ce or biomass of differenpan>t decomposers is of critical importanpan>ce for decomposition (e.g. Bergfur et al., 2007; Ebelinpan>g et al., 2014; Manpan>ninpan>g anpan>d Cutler, 2018). Evenpan> at constanpan>t richness anpan>d communpan>ity composition, strong decreases inpan> abunpan>danpan>ce canpan> have importanpan>t impacts on ecosystem funpan>ctioninpan>g (Spaak et al., 2017; but see Dainpan>ese et al., 2019). It is beyond the scope of the presenpan>t meta-anpan>alysis to disenpan>tanpan>gle the effects of biodiversity from the effects of abunpan>danpan>ce, anpan>d we founpan>d that both contributed to explainpan> shifts inpan> decomposition inpan> separate anpan>alyses. It is inpan>terestinpan>g to note that the few cases where negative effect sizes of pan> class="Disease">stressors on biodiversity were associated with positive effect sizes on decomposition were also cases where decomposer abundance was positively associated with stressors (Figure 4). Although we cannot specifically test this with the present data, it seems that in those particular cases (Lucisine et al., 2015), increases in decomposer abundance counteracted the negative effects of decreases in decomposer diversity (Dornelas et al., 2019). Those results could therefore be in line with the mass-ratio hypothesis (Grime, 1998; Smith and Knapp, 2003). Indeed, an exclusion experiment showed that dominant, small, detritivores can compensate reductions in litter decomposition caused by the removal of large detritivores (Cárdenas et al., 2017). These concomitant shifts in both diversity and abundance further have important implications for our estimates of diversity responses, as studies mostly reported richness to estimate decomposer diversity, but rarely corrected for the sampling effort (Gotelli and Colwell, 2001). This means that lower abundances rather than a lower number of species per se might have directly caused some of the negative effects on biodiversity reported here (Chase and Knight, 2013). This common caveat in meta-analysis approaches that rely on how individual studies report biodiversity, also applies to the present study, and reinforces the importance of reporting raw data in future studies on the impacts of chemical stressors on biodiversity and ecosystem functioning.

The effects of changes in decomposer diversity and abundance on decomposition found in the present study might also have channeled changes in community and food-web structure not captured by our biodiversity metrics. Changes in keystone species (Hättenschwiler et al., 2005), functional diversity (Cadotte et al., 2011; Dangles et al., 2012; Heemsbergen et al., 2004), vertical diversity (Gessner et al., 2010; Melguizo-Ruiz et al., 2020; Wang and Brose, 2018; Zhao et al., 2019), or dominance patterns (Dangles and Malmqvist, 2004) might have shifted concomitantly to taxonomic diversity and abundance. Moreover, these different components of diversity might act at different timings of decomposition (Oliveira et al., 2019). Unfortunately, studies rarely reported such measurements together with decomposition. For example in our dataset, only seven studies reported evenness. Future studies need to explore shifts in decomposer community composition in more detail to better understand what particular aspect of biodiversity is responsible for changes in decomposition rates (Giling et al., 2019; Hättenschwiler et al., 2005). In particular, few of the included studies reported comparable functional groups allowing to address the effect of functional diversity across the multiple systems and taxonomic groups addressed by the present analysis. Future synthesis work could specifically address the effect of functional diversity, by focusing on a given system type. Indeed, there is ample evidence that shifts in functional diversity are crucial for decomposition (Heemsbergen et al., 2004), and that facilitative interactions occur primarily between decomposers of contrasting body size (Dangles et al., 2012; Tonin et al., 2018). This is especially the case for interactions between animal and microbial decomposers, where fragmentation of litter by detritivores facilitates access for microbial decomposers (Eisenhauer et al., 2010; Hättenschwiler et al., 2005; Yang et al., 2012).

Here, we found that invertebrates were more affected by chemical stressors thanpan> microbes, across aquatic anpan>d terrestrial ecosystems. Invertebrate decomposers are particularly senpan>sitive to the impacts of pan> class="Chemical">metals and pesticides (Hogsden and Harding, 2012; Pelosi et al., 2014; Peters et al., 2013; Schäfer, 2019). Microbial decomposers are known to be sensitive to metals (Giller et al., 2009) and pesticides as well (DeLorenzo et al., 2001). Nevertheless, our result is consistent with the general expectation that larger organisms are more sensitive to environmental change due to longer generation time, higher energetic demands and lower population densities (Hines et al., 2015; Sheridan and Bickford, 2011; Woodward et al., 2005; Yvon-Durocher et al., 2011; Baas and Kooijman, 2015). These different sensitivities between groups of decomposers could imply that the biodiversity-mediated effects of stressors on decomposition are more strongly linked to shifts in invertebrates than microbes, as reported in a previous review (Peters et al., 2013). However, in another meta-analysis focusing on microbial-driven decomposition rates, changes in fungal biomass and richness explained shifts in decomposition under the impacts of chemical stressors, but also of nutrient enrichment (Lecerf and Chauvet, 2008).

Nutrient-induced changes in decomposition were not related to shifts in decomposer diversity

The impacts of n class="Gene">nutrienpan>t enpan>richmenpan>t on litter decomposition anpan>d decomposer diversity were differenpan>t from those caused by pan> class="Disease">stressors, confirming our expectations. These different biodiversity and function responses led to different emergent relationships between decomposer diversity and decomposition compared to stressors. We found that nutrients had a variety of effects ranging from positive to negative depending on the taxonomic group (Figure 7) and nutrient intensity (Figure 6), and resulting in neutral overall mean effects (Figure 3). Previous syntheses also found positive (Ferreira et al., 2015) as well as inconsistent (Knorr et al., 2005) responses of decomposition rates to nutrient enrichment in streams. The relatively small mean effect of nutrient enrichment on decomposition in the present meta-analysis could be explained by the use of correlation as an effect size, which does not capture potentially non-monotonic responses of decomposition to nutrients (Woodward et al., 2012). However, we noted that most of the studies included in the present meta-analysis did not individually span nutrient gradients sufficiently large to capture this potential non-monotonous response. Taken together, the studies show positive effects on decomposition at low nutrient intensities that shifted toward neutral to negative effects at higher intensities (Figure 6), which is consistent with previous findings (Ferreira et al., 2015; Woodward et al., 2012). Low-nutrient intensities might have enhanced microbial activity and biomass by alleviating resource limitation, resulting in enhanced decomposition. At higher intensities, however, negative impacts on invertebrates might have decreased decomposition rates (Peters et al., 2013; Woodward et al., 2012).

These n class="Gene">nutrienpan>t inpan>tenpan>sity patternpan>s contrasted with the results for chemical pan> class="Disease">stressors. The overall negative effects of stressors (Figure 3) on decomposition were not explained by stressor intensity levels (Figure 6), and there was mixed support for a stressor intensity effect on decomposer diversity based on two complementary data analysis approaches (SEM based on data resampling (Figure 5) vs. second level meta-analysis Figure 6). Thus, negative responses to chemical stressors happened across the range of stressor intensity. Such contrasting patterns between stressor and nutrient intensity effects may reflect the greater number of stressor types (different metals, pesticides, mixtures) covered by individual studies compared to the limited number of nutrients. In addition, due to the higher variability of stressor types, we relied on more variable sources to standardize stressor levels compared to nutrients in the diversity dataset (Materials and methods, Appendix 1—table 1). With the data at hand, it was not possible to test the influence of the environmental quality criteria used to standardize stressor and nutrient levels, because such an effect would be confounded with stressor or nutrient types. The datasets were all dominated by environmental quality criteria based on similar methodologies (for 75% to 100% of observations, see Material and Methods). However, future studies focusing on stressor intensity effects across ecosystems would greatly benefit from coordinated efforts to derive quality criteria encompassing the vast and rapidly increasing number of chemical stressors (Wang et al., 2020).

Appendix 1—table 1.

Environmental quality criteria for stressors and nutrients.

Quality criteria were used to standardized the intensity levels of the different chemical stressors across studies included in the meta-analysis.

System	Chemical or nutrient	Unit1	Unit2	Quality criteria	Citation
aquatic	fungicide: pyrimethanil	µg/l	-	0.69	Abelho M, Martins TF, Shinn C, Moreira-Santos M, Ribeiro R. 2016. Effects of the fungicide pyrimethanil on biofilm and organic matter processing in outdoor lentic mesocosms. Ecotoxicology 25:121–131.
aquatic	fungicide: tebuconazole	µg/l	-	0.10	https://echa.europa.eu/documents/10162/41e9d7aa-4559-f904-9cb5-0a0d5f0d6445
aquatic	As	µg/l	-	13.00	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.316
aquatic	Al	µg/l	-	87.00	https://www.govinfo.gov/content/pkg/FR-2018-12-21/pdf/2018-27745.pdf
aquatic	Cu	µg/l	-	10.10	https://echa.europa.eu/brief-profile/-/briefprofile/100.124.825
aquatic	Zn	µg/l	-	20.60	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.341
aquatic	Fe	µg/l	-	1000.00	https://www.epa.gov/wqc/national-recommended-water-quality-criteria-aquatic-life-criteria-table
aquatic	Mn	µg/l	-	1000.00	https://www.epa.gov/wqc/national-recommended-water-quality-criteria-aquatic-life-criteria-table
aquatic	Hg	µg/l	-	0.06	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.278
aquatic	Cd	µg/l	-	0.19	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.320
aquatic	insecticide: chlorpyrifos	µg/l	-	0.08	https://www.epa.gov/wqc/national-recommended-water-quality-criteria-aquatic-life-criteria-table
aquatic	phenanthrene	µg/l	-	51.40	Wu, J. Y., Yan, Z. G., Liu, Z. T., Liu, J. D., Liang, F., Wang, X. N., & Wang, W. L. (2015). Development of water quality criteria for phenanthrene and comparison of the sensitivity between native and non-native species. Environmental Pollution, 196, 141-146.
aquatic	Zn	mg/kg	-	117.80	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.341
aquatic	Cd	mg/kg	-	1.80	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.320
aquatic	Hg	mg/kg	-	9.30	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.278
aquatic	Pb	mg/kg	-	186.00	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.273
terrestrial	Cu	mg/kg	-	106.35	https://echa.europa.eu/brief-profile/-/briefprofile/100.124.825
terrestrial	Zn	mg/kg	-	35.60	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.341
terrestrial	Ni	mg/kg	-	29.90	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.283
terrestrial	Mn	mg/kg	-	3.40	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.277
terrestrial	Hg	µg/kg	-	22.00	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.278
terrestrial	Pb	mg/kg	-	212.00	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.273
terrestrial	Cd	mg/kg	-	0.90	https://echa.europa.eu/brief-profile/-/briefprofile/100.028.320
terrestrial	insecticide: chlorpyrifos	kg/ha	-	1.25	Iwai CB, Noller B. 2010. Ecotoxicological assessment of diffuse pollution using biomonitoring tool for sustainable land use in Thailand. Journal of Environmental Sciences 22:858–863.
terrestrial	insecticide: endosulfan	kg/ha	-	1.25	Iwai CB, Noller B. 2010. Ecotoxicological assessment of diffuse pollution using biomonitoring tool for sustainable land use in Thailand. Journal of Environmental Sciences 22:858–863.
terrestrial	herbicide: atrazine	kg/ha	-	1.88	Iwai CB, Noller B. 2010. Ecotoxicological assessment of diffuse pollution using biomonitoring tool for sustainable land use in Thailand. Journal of Environmental Sciences 22:858–863.
terrestrial	insecticide: carbofuran	kg/ha	-	31.25	Iwai CB, Noller B. 2010. Ecotoxicological assessment of diffuse pollution using biomonitoring tool for sustainable land use in Thailand. Journal of Environmental Sciences 22:858–863.
aquatic	pesticide mixture	arbitrary	-	1.00	Talk A. 2016. Effects of multiple butLow pesticide loads on aquatic fungal communities colonizing leaf litter. Journal of EnvironmentalSciences 46:116–125.
terrestrial	herbicide: glyphosate	kg/ha	-	4.32	European Food Safety Authority (EFSA). Conclusion on the peer review of the pesticide risk assessment of the active substance glyphosate. EFSA Journal 13, (2015).
terrestrial	herbicide: simazine	kg/ha	-	0.10	https://ec.europa.eu/food/plant/pesticides/eu-pesticides-database/public/?event = activesubstance.detail and language = EN and selectedID = 1853
aquatic	pesticide mixture	sum or max of TU (toxic units)	-	−3.50	Schäfer, R. B., Caquet, T., Siimes, K., Mueller, R., Lagadic, L., & Liess, M. (2007). Effects of pesticides on community structure and ecosystem functions in agricultural streams of three biogeographical regions in Europe. Science of the Total Environment, 382(2-3), 272-285.
aquatic	DIN	mg/l	N	3.05	Ministère de l’Environnement, de l’Énergie et de la Mer. Guide technique Relatif à l’évaluation de l’état des eaux de surface continen- tales (cours d’eau, canaux, plans d’eau). (2016).
aquatic	NH4+	mg/l	NH4	0.10	Ministère de l’Environnement, de l’Énergie et de la Mer. Guide technique Relatif à l’évaluation de l’état des eaux de surface continen- tales (cours d’eau, canaux, plans d’eau). (2016).
aquatic	NO3	mg/l	NO3	10.00	Ministère de l’Environnement, de l’Énergie et de la Mer. Guide technique Relatif à l’évaluation de l’état des eaux de surface continen- tales (cours d’eau, canaux, plans d’eau). (2016).
aquatic	NO2	mg/l	NO2	0.10	Ministère de l’Environnement, de l’Énergie et de la Mer. Guide technique Relatif à l’évaluation de l’état des eaux de surface continen- tales (cours d’eau, canaux, plans d’eau). (2016).
aquatic	Total_N	mg/l	N	0.67	US EPA, O. Water Quality Criteria. US EPA (2013). Available at: https://www.epa.gov/wqc. (Accessed: 7th January 2019)
aquatic	SRP	mg/l	PO43	0.10	Guide technique Relatif à l’évaluation de l’état des eaux de surface continen- tales (cours d’eau, canaux, plans d’eau). (Ministère de l’Environnement, de l’Énergie et de la Mer, 2016).
aquatic	Total_P	mg/l	P	0.05	Guide technique Relatif à l’évaluation de l’état des eaux de surface continen- tales (cours d’eau, canaux, plans d’eau). (Ministère de l’Environnement, de l’Énergie et de la Mer, 2016).
terrestrial	N deposition	kg/ha/yr	N	20.00	Pardo, L.H., Fenn, M.E., Goodale, C.L., Geiser, L.H., Driscoll, C.T., Allen, E.B., Baron, J.S., Bobbink, R., Bowman, W.D., Clark, C.M., Emmett, B., Gilliam, F.S., Greaver, T.L., Hall, S.J., Lilleskov, E.A., Liu, L., Lynch, J.A., Nadelhoffer, K.J., Perakis, S.S., Robin-Abbott, M.J., Stoddard, J.L., Weathers, K.C. and Dennis, R.L. (2011), Effects of nitrogen deposition and empirical nitrogen critical loads for ecoregions of the United States. Ecological Applications, 21: 3049-3082. doi:10.1890/10-2341.1; derived critical loads (i.e. level of deposition below which no detrimental ecological effect occurs over the long term according to current knowledge) from empirical data for various (plant) species and ecosystems
terrestrial	P fertilization	kg/ha/yr	P	35.00	Amery, F., & Schoumans, O. F. (2014). Agricultural phosphorus legislation in Europe. Institute for Agricultural and Fisheries Research (ILVO).

Contrary to our expectation, n class="Gene">nutrienpan>t-inpan>duced shifts inpan> decomposer diversity anpan>d abunpan>danpan>ce were not associated with shifts inpan> decomposition rates across studies. We founpan>d that inpan>creasinpan>g pan> class="Gene">nutrient intensity decreased the effects on decomposition and on decomposer diversity, but not on decomposer abundance. Statistically controlling for the effect of nutrient intensity with SEM indicated no residual association between shifts in decomposer diversity or abundance and in decomposition rates, that is a non-significant BEF relationship. Changes in microbial abundance in response to nitrogen deposition explained the responses of different ecosystem functions in terrestrial systems in previous meta-analyses (García-Palacios et al., 2015; Treseder, 2008). Here, we show that this pattern cannot be generalized across aquatic and terrestrial systems and across animal and microbial decomposers. Contrary to stressors, when the diversity and abundance of animal and microbial decomposers were not affected by nutrients, we observed large positive and negative shifts in decomposition (intercepts of Figure 4), that were explained by nutrient intensity (Figure 4: negative effects on decomposition at invariant biodiversity are associated with high intensities and positive effects with lower intensities). Together, these results show that nutrient-induced shifts in decomposer diversity were not as strong drivers of decomposition changes as stressor-induced biodiversity shifts. These differences may be partly due to the different mechanisms underlying the effects of stressors and nutrients. Based on previous studies, we speculate that our results are due to the complex responses of animal and microbial decomposers at different nutrient intensities (Ferreira et al., 2015; Lecerf and Chauvet, 2008; Treseder, 2008; Woodward et al., 2012).

Animal decomposers showed a stronger response to n class="Gene">nutrienpan>ts thanpan> microbes. Invertebrate decomposers overall decreased inpan> diversity, but they inpan>creased inpan> abunpan>danpan>ce unpan>der pan> class="Gene">nutrient enrichment. These results could reflect a loss of sensitive taxa to the benefit of tolerant taxa that were able to use additional resources and would then increase in density (Bergfur et al., 2007). Overall, microbial decomposers responded little to nutrient enrichment, probably reflecting a mixture of positive and negative effects that nutrients can have on microbial growth (Lecerf and Chauvet, 2008; Treseder, 2008), as well as on different microbial taxa. Indeed, nutrients can alleviate resource limitations at low intensities, but can also exert toxic effects at high intensities. The initial levels of nutrients thus condition subsequent responses in decomposers and decomposition to nutrient enrichment (Ferreira et al., 2015; Knorr et al., 2005). Furthermore, at high intensities, nutrients can be associated with other chemical stressors (e.g. pesticides in agricultural runoffs) (Ferreira et al., 2015; Woodward et al., 2012). The influence of interactive effects of stressors and nutrients was impossible to quantify with the data at hand, given that only a few experiments assessed the effects of both drivers independently, but many observational studies may have been confounded by such joint effects. Chemical stressors and nutrients are often co-occurring in e.g. agricultural landscapes, and the consequences of such combinations are still poorly understood (Alexander et al., 2013; Alexander et al., 2016; Barmentlo et al., 2018; Chará-Serna et al., 2019; Chará-Serna and Richardson, 2018; Fernández et al., 2016). Furthermore, stressor and nutrient effects might be modulated by climatic and other environmental conditions, and studies on interaction effects are scarce (Rillig et al., 2019; Thakur et al., 2018). Finally, although our comparison of stressors versus resources allowed us to test a clear concept, any kind of grouping in ecological studies may mask some of the variation within the categories and future studies may be interested in different categories. Indeed, a given environmental change driver can represent a stressor for a given species, and a resource for another species (Connell et al., 2018). As data availability improves, future work could include different environmental change drivers. This would also allow to test additional groupings of drivers and ecological concepts unifying stressors and resources (De Laender, 2018; Harley et al., 2017).

Conclusions

This study brings new insights into the real-world patterns relating ecosystem function to non-random changes in biodiversity induced by environmental change. We found that the consequences of changes in biodiversity for ecosystem functioning depend on the type of environmental change. Real-world scenarios do not necessarily involve concomitant changes in both biodiversity and function across terrestrial and aquatic systems. We further found that with the environmental quality criteria used in risk assessment, there were already significant positive and negative effects on decomposers and decomposition (Figure 6), highlighting the need to better incorporate biodiversity and ecosystem function into ecological risk assessment programs (De Laender and Janssen, 2013). Finally, we report overall negative effects of chemical stressors on biodiversity anpan>d ecosystem funpan>ctioninpan>g across terrestrial anpan>d aquatic ecosystems that reinpan>force recenpan>t calls to consider chemical pan> class="Disease">stressors as important global change drivers and address their impacts on biodiversity and ecosystems (Bernhardt et al., 2017; Mazor et al., 2018; Steffen et al., 2015). Positive real-world BEF relationships may be particularly significant in cases where environmental changes decrease biodiversity, such as in the case of chemical stressors. Such information are crucial if we are to design policy and conservation strategies able to reconcile human development with biodiversity conservation.

Materials and methods

Data collection

We searched the Web of Science for studies that addressed the impact of environmental drivers and recorded decomposer community responses and litter decomposition rates. The search strategy is fully reported in Supplementary Methods (Appendix 1). The search retrieved 2536 references. Abstracts and titles were screened to identify a final set of 61 records that met our inclusion criteria (PRISMA plot, Appendix 1—figure 1, and list of included references (Appendix 4). To be included in the meta-analysis, studies had to:

Appendix 1—figure 1.

PRISMA plot describing the data collection steps of the meta-analysis.

SEM = structural equation modeling.

Report litter decomposition (rates, mass loss, proportion of mass remaining) and the diversity, abundance, or biomass of decomposers at sites differing in chemical n class="Disease">stressor or n class="Gene">nutrient levels.

Focus on naturally assembled communities subjected to the impact of chemical n class="Disease">stressors or pan> class="Gene">nutrient enrichment. Studies that manipulated decomposer diversity directly were not considered to only focus on non-random biodiversity change scenarios. We included mesocosm studies only when they used field-sampled communities and left time for the community to reach an equilibrium in mesocosms in order to reflect real-world conditions as much as possible.

Report the response of animal (benthic macroinvertebrates, or soil micro, meso or macrofauna) or microbial decomposers (bacteria or fungi from decomposing leaves or in surrounding n class="Chemical">water or soil samples).

Report decomposer abundance (density or biomass), or decomposer diversity (taxa richness, Shannon diversity, evenness).

When a reference reported different environmental change drivers or geographical areas with a specific reference site for each case, we considered these as individual (case) studies (García-Palacios et al., 2015). We extracted means or sums, standard deviations, and sample sizes of litter decomposition, decomposer diversity, and abundance (outcomes) in non-impacted vs. impacted sites (control-treatment studies), or at each site when gradients of chemical stressors or pan> class="Gene">nutrients were investigated (gradient studies). When response variables were reported at different time points, we kept only the last time point to capture long-term responses. For studies reporting decomposition, decomposer abundance or diversity for several litter types (e.g. different litter species), several groups of organisms (e.g. functional feeding groups for macroinvertebrates), and several diversity metrics (e.g. Shannon indices and taxon richness), we created separate observations within case studies. We also extracted chemical stressor or nutrient levels at those sites (water, soil, or sediment concentrations of chemical stressors or nutrients, or application rate of pesticides or fertilizers). The study type (experimental vs. observational), taxonomic group (animal decomposers or microbial decomposers) and metric of diversity (taxa richness or diversity indices (Shannon diversity and evenness)) were also recorded. We used the online software Webplotdigitizer to extract data from figures (Rohatgi, 2018). We converted standard errors and confidence intervals into standard deviations using the equations in Lajeunesse, 2013. When reported as mass loss, litter decomposition data were transformed into k rates using the exponential decay equation used in Ferreira et al., 2015.

Effect size calculation

We used z-transformed correlation coefficients as effect sizes in order to cope with the heterogeneity of data and study types (Koricheva et al., 2013). For control-treatment studies, we first calculated Hedge’s d, and then transformed Hedge’s d into correlation coefficients (Lajeunesse, 2013). For gradient studies (four or more treatment levels), we calculated correlation coefficients between the mean values of abundance, diversity, or decomposition rate and the corresponding chemical n class="Disease">stressor or pan> class="Gene">nutrient concentrations. When means, standard deviations, or sample sizes were missing, we contacted the authors to retrieve the data. When the information could not be retrieved, standard deviations were approximated from the data, using the linear relationship between mean values and standard deviations across our datasets (Lajeunesse, 2013).

Standardization of chemical stressors and nutrient enrichment intensities

Given the variability in the different stressors anpan>d pan> class="Gene">nutrients combinations in the studies, stressor and nutrient levels were standardized into a common environmental change driver intensity () as follows:where were environmental quality criteria set by European or US environmental authorities for the chemical stressor or nutrient considered (Appendix 1—table 1), and were the concentrations of the chemical stressor or nutrient at the treatment or impacted sites. When multiple stressors or nutrients were reported, we used the standardized intensity of the stressor or nutrient corresponding to the highest standardized intensity for the rest of the analyses.

We used consistent sources for the environmental quality criteria as much as possible. For chemicals, we relied primarily on quality criteria from the European Chemical Agency (ECHA) and United States Environmental Protection Agency (USEPA) that use standardized procedures across aquatic and terrestrial realms based on ecotoxicological data. For n class="Gene">nutrienpan>ts, we relied mostly on Europeanpan> pan> class="Chemical">Water Framework Directive (WFD) benchmarks. Using various sources for those quality criteria was inevitable due to the high number of chemicals and the various way the authors reported stressor or nutrient levels in individual studies. When we could not find quality criteria for the stressors or nutrients considered in the studies in our main sources, we relied on the authors’ statements and expert knowledge regarding their stressor or nutrient levels (e.g. citation for ecotoxicological data, or synthesis studies, or recommended application rates of pesticides [Appendix 1—table 1]). Despite this, the final datasets were all dominated by similar sources for standardizing stressor and nutrient intensity levels: thresholds from ECHA or USEPA for 80% and 90% of observations in the stressor-diversity and stressor-abundance datasets, respectively, and for nutrients, thresholds from WFD for 100% and 75% of observations in the nutrient-diversity and nutrient-abundance datasets, respectively.

Overall effects of chemical stressors and nutrient enrichment: first-level meta-analysis

We first tested the differences between the effects of chemical stressors anpan>d pan> class="Gene">nutrient enrichment on decomposer diversity, abundance and litter decomposition responses by quantifying the grand mean effect sizes on the three response variables (first level meta-analysis). Three separate meta-analyses were conducted, one for each response variable, and included the type of driver (stressors or nutrients) as a categorical moderator, and a random effect of the case study. We used a weighted meta-analysis giving more weight to effect sizes derived from studies with larger sample sizes. Weights were the inverse of the variance in z-transformed correlation coefficients (Viechtbauer, 2010). Publication bias was evaluated using funnel plots with environmental change driver type as covariate. The intercepts from Egger’s regressions (standardized effect size vs. precision = 1/SE) were inspected for significant deviation from zero that would indicate publication bias (Koricheva et al., 2013). Residual plots were used to detect strong deviation from normality and outliers. We estimated the grand mean effect sizes and compared the effect of chemical stressors and of nutrients using Wald-type chi-square tests. The rma.mv() function of the R package metafor was used (R Development Core Team, 2018; Viechtbauer, 2010).

Relationship between biodiversity and decomposition: Structural equation modelling

An SEM was fitted to estimate the relationship between decomposer diversity or abundance and litter decomposition responses to environmental change drivers while controlling for the joint influence of n class="Disease">stressor or pan> class="Gene">nutrient intensity and categorical covariates. We used piecewise SEM (Lefcheck, 2001) estimating two linear mixed effect models, one for decomposition () and one for decomposer diversity or abundance responses (), with a random effect of the case study on the intercepts. These two sub-models embedded in the piecewise SEM were the second-level meta-analyses in our hierarchical approach. The random effect structure, weighting approach and variance structure were coded with the R package nlme (Pinheiro et al., 2018) in a way that fully reproduced the meta-analysis approach of weighting and of known residual variance (Viechtbauer, 2017):

This SEM was tested separately for each of four datasets: n class="Disease">Stressors – Biodiversity; pan> class="Disease">Stressors – Abundance; Nutrients – Biodiversity and Nutrients – Abundance datasets. The influence of the diversity metric (diversity indices versus taxa richness) was tested in the Biodiversity datasets only. We initially considered more complex model structures, but were unable to use them for analysis due to data limitations (in particular the effect of the ecosystem type and of interactions between our covariates).

Outliers, relationships between covariates, and non-linear patterns between continuous covariates were explored graphically. Studies often reported different decomposer diversity or abundance values for the same litter decomposition (e.g. when several taxonomic or functional groups were reported in the same litterbag). This variability could have affected the model estimates. We thus used data resampling to account for duplicated effect sizes on litter decomposition in the analyses. A stratified resampling was conducted, where for each duplicated value of effect size on decomposition, one randomly selected effect size on biodiversity was kept at each out of 1000 iterations. The models were fitted for each data resampling iteration, and we averaged model estimates and statistics across iterations and used the means as final values (path coefficients and standard error of the path and intercepts, Chi-square statistics and AICs).

Goodness-of-fit of the SEMs was assessed using directed separation tests based on the Fisher’s C statistic. We used mediation tests to explore the significance of the path between decomposer diversity or abundance and litter decomposition based on the Fisher’s C statistic of SEM that did not include the biodiversity-mediated path (Lefcheck, 2001; Shipley, 2009). We calculated the p-value associated with the mean Fisher’s C statistic across data resampling iterations (p-value<0.05 indicated poor model fit). The AICs of models with and without the biodiversity-mediated paths were further compared using averaged AICs across data resampling iterations. We considered the biodiversity (or abundance) path to be consistent with the data when the SEM without the biodiversity-path had p-value<0.05 (poor fit) and was not associated with a better AIC value (i.e. lower than two units) than the SEM including the biodiversity path. Residuals from the two sub-models of each SEM were graphically evaluated for strong departure to normality and relationship with the fitted values (Duffy et al., 2015). For these analyses, we averaged the residuals across data resampling iterations for each observation. We finally compared the relative magnitude of the biodiversity-mediated path versus the direct path from n class="Disease">stressor or pan> class="Gene">nutrient intensity to litter decomposition based on the mathematical product of the standardized path coefficients (Grace, 2006).

Moderator analyses: second-level meta-analyses

In order to quantify the influence of the categorical (study type, taxonomic group and diversity metrics) and continuous (environmental change intensity) moderators on the three response variables, we further analyzed the results of the second-level meta-analyses (i.e. the sub-models embedded in the SEMs). The data resampling used in the SEM was no longer necessary, because there were no repeated values of decomposition matching different decomposer diversity or abundance measurements in this univariate approach. We quantified the effects of the different moderators based on the Wald-type chi-square tests derived with the R package metafor (Viechtbauer, 2010).

Sensitivity analyses

We finally tested the robustness of the results to the approximation of standard deviations, the presence of extreme values, and the metric of effect size used. The analyses were re-run with datasets that did not include the effect sizes for which we approximated standard deviations, for datasets that did not include extreme values of effect sizes (values beyond the whiskers of boxplots that is below quantile 1 minus 1.5 times the interquartile range or above quantile 3 plus 1.5 times the interquartile range). Finally, we calculated log-response ratios instead of correlation coefficients as effect sizes and re-run the analyses.

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Many ecosystems are severely affected by n class="Species">human activity, yet we still knpan>ow surprisinpan>gly little about the modulatinpan>g role of specific elemenpan>ts on ecosystem health on a global scale. This paper provides a convinpan>cinpan>g anpan>d timely synpan>thesis of the effects of chemical pan> class="Disease">stressors and nutrient enrichment on soil biodiversity, and a major ecosystem process they modulate – litter decomposition. This work will inspire new research linking ecological stoichiometry to ecosystem services.

Decision letter after peer review:

Thank you for submitting your article "Biodiversity mediates the effects of stressors but not nutrients on litter decomposition" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Christian Rutz as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reports with one another, and the Reviewing Editor has drafted this decision letter to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to n class="Disease">COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are askinpan>g editors to accept without delay manpan>uscripts, like yours, that they judge canpan> stanpan>d as eLife papers without additional data, even if they feel that they would make the manpan>uscript stronger. Thus the revisions requested below only address clarity anpan>d presentation.

Summary:

The paper describes a meta-analysis on the effects of n class="Gene">nutrienpan>ts anpan>d chemical pan> class="Disease">stressors on biodiversity and the relationship with ecosystem functioning. Using two meta-analyses and Structural Equation Modelling, they find that chemical stressors (but not added nutrients) affect leaf litter decomposition, through changes in biodiversity and abundance of organisms. Interestingly, the authors do not establish a stressor intensity-effect relationship. They find minor effects of nutrients on the biodiversity-ecosystem functioning relationship, though they establish relationships with nutrient intensity, depicting real-world scenarios of human impacts on ecosystems.

Essential revisions:

1) The authors run different analyses – for the average reader it may be difficult to understand why there is a stressor effect inpan> the SEM, but no pan> class="Disease">stressor intensity-effect relationship. We find this aspect under-appreciated in the Discussion, and it could be described more clearly what the reasons for this difference are. In addition, it should be discussed to which extent the variable basis for the values used to characterize stressor intensity may be the reason for this finding. There is a relationship for nutrients, where the authors have a relatively consistent data source, whereas, for chemical stressors where the data sources vary, including the methods used to derive quality criteria, they found no relationship. Our suspicion is therefore that this is partly due to methodical reasons. To clarify, problems may arise if you choose a benchmark Y for stressor S1, but multiple benchmarks Zi where i from 1 to n for Stressor S2i that are not consistent with each other. Now imagine you find a relationship for one but not for the other with stressor intensity – this can simply be because you have variable benchmarks Zi and is not necessarily due to the stressor.

2) The stressor definpan>ition should be reconsidered. The authors definpan>e: "pan> class="Disease">stressors (e.g., temperature, drought, chemicals) and resource enrichment (e.g., of CO2 or mineral nutrients)." From a theoretical perspective, the temperature would also be a niche-defining resource. From a practical perspective, it is confusing for the reader to contrast "chemical stressors" to CO2 and NO3, the latter of which are clearly chemicals and often anthropogenic. Then both are "environmental change drivers", which makes the terminology even more complicated for the reader. Why not just use stressors and environmental change drivers, which have different niche-defining properties? Some are unimodal (nutrients, temperature), others exhibit a threshold relationship and others may be log-linear.

Essential revisions:

1) The authors run different analyses – for the average reader it may be difficult to understand why there is a stressor effect inpan> the SEM, but no pan> class="Disease">stressor intensity-effect relationship. We find this aspect under-appreciated in the Discussion, and it could be described more clearly what the reasons for this difference are.

We agree that the Discussion was not clear enough regarding this important point. We conducted two approaches for analyzing the data: an SEM and a meta-regression. The results are slightly different because different datasets were used. For the SEM, we used data resampling. We tested 1000 iterations of the model, each iteration based on a random sample of the effect sizes on biodiversity or abundance. This approach was done to account for the fact that studies often reported multiple measures of biodiversity or abundance response in the same litterbag where litter decomposition was assessed. Thus, the data had many cases where a unique effect size on decomposition was associated with multiple effect sizes on biodiversity or abundance. Such a data structure could have biased the coefficients derived from the SEM, and we conducted data resampling to tackle this issue.

In the meta-regression depicted in Figure 6, we show the correlations between biodiversity responses and n class="Disease">stressor inpan>tensity derived from our second-level meta-anpan>alysis. This is not accountinpan>g for the data structure (all the effect sizes on biodiversity or abundanpan>ce are represented anpan>d inpan>cluded inpan> the regression). The data resamplinpan>g used inpan> the SEM was no longer necessary because there were no repeated values of decomposition matchinpan>g different decomposer diversity or abundanpan>ce measurements inpan> this univariate approach. This approach further complemented the inpan>formation provided by the SEM by testinpan>g for the effect of moderators (such as study type anpan>d taxonomic groups). Therefore, the results from SEM anpan>d meta-regression differ because different datasets are used.

Furthermore, it is important to note that the magnitude of the relationships between Figure 5 (SEM) and Figure 6 (n class="Disease">stressor inpan>tenpan>sity effect) are not comparable. Figure 5 shows the stanpan>dardized coefficienpan>ts, while Figure 6 shows unpan>stanpan>dardized coefficienpan>ts. The unpan>stanpan>dardized coefficienpan>ts from the SEM (Appenpan>dix 2—table 2) anpan>d the unpan>stanpan>dardized coefficienpan>ts of the pan> class="Disease">stressor intensity relationships (Figure 6) are actually not that different from each other (with the SEM the slope was -0.10 (SE: 0.04, p = 0.02), while the meta-regression slope was = -0.05 (SE: 0.03, p = 0.11)). Again, these differences can be explained by the fact that different data underlie those results.

Finally, we re-tested the meta-regression between chemical n class="Disease">stressor inpan>tenpan>sity anpan>d decomposer diversity responses, this time inpan>cludinpan>g anpan> additional ranpan>dom effect of the litterbag (therefore addressinpan>g the data structure similarly as inpan> the SEM). This approach yielded more similar results as the SEM (a signpan>ificanpan>t pan> class="Disease">stressor intensity effect on diversity responses: slope = -0.09, p=.009). We initially tested if including such an additional random effect improved the models. As this was not the case (Likelihood ratio test: 2.05, p = 0.15, ΔAIC = 0.05), we decided to keep the random effect structure minimal, thereby not including this additional litterbag random effect. The litterbag effect might have been confounded by the diversity metric effect in our models given that studies often reported different diversity metrics for the same litterbag sample.

We modified the main text to describe the potential reasons for this difference with more details:

“There was mixed support for a n class="Disease">stressor inpan>tenpan>sity effect on decomposer diversity across the two approaches: decomposer diversity responses decreased with pan> class="Disease">stressor intensity according to the SEM (Figure 5) but this trend was not significant according to the second level meta-analysis (Figure 6). […] The differences between the two approaches can be explained by the different data involved.”

We also discuss these findings in the Discussion section:

“These nutrient inpan>tensity patternpan>s contrasted with the results for chemical pan> class="Disease">stressors. […] Thus, negative responses to chemical stressors happened across the range of stressor intensity.”

In addition, it should be discussed to which extent the variable basis for the values used to characterize stressor inpan>tenpan>sity may be the reason for this finpan>dinpan>g. There is a relationship for pan> class="Gene">nutrients, where the authors have a relatively consistent data source, whereas, for chemical stressors where the data sources vary, including the methods used to derive quality criteria, they found no relationship. Our suspicion is therefore that this is partly due to methodical reasons. To clarify, problems may arise if you choose a benchmark Y for stressor S1, but multiple benchmarks Zi where i from 1 to n for Stressor S2i that are not consistent with each other. Now imagine you find a relationship for one but not for the other with stressor intensity – this can simply be because you have variable benchmarks Zi and is not necessarily due to the stressor.

This is a very good point that we carefully considered prior to our analyses. We fully agree that discussing the variable standardization used is important, and we modified the text accordingly (see below). Chemical stressors inpan> our meta-anpan>alysis inpan>clude a wide ranpan>ge of chemicals (differenpan>t pan> class="Chemical">metals, and pesticides). On the opposite, nutrient additions were far less variable, with only a hand full of nutrient types covered (NH4, NO3, NO2, PO43-, total N and total P). It would be interesting to test the robustness of the results if we had only incorporated observations for which similar threshold sources could be derived. However, such a sensitivity analysis would not separate the “variable threshold effect” from the effect of certain type of chemicals versus others. This is because variable thresholds had to be used for particular types of pollutants (pesticide use, pesticide mixtures see below). Therefore, we are unsure whether the methodology to standardize stressor and nutrient intensity is the reason for the different relationships that we found. We rather think that the variable nature of chemical stressors versus nutrient enrichment, as defined in our analysis, is the reason for those differences.

The quality criteria values originated from consistent sources as much as possible. For chemicals, we mostly relied on two main authoritative sources: the European chemicals agency (ECHA) and the United States Environmental Protection Agency (USEPA). In our final datasets, we were able to retrieve quality criteria from those standard sources for 80% and 90% of the observations for diversity and abundance datasets, respectively. We preferably used ECHA and USEPA, because they derive quality criteria in a standardized way. Both are using similar approaches that rely on ecotoxicological data compiled for as many taxa as possible for different types of ecosystems. Importantly, these methodologies are consistent across aquatic and terrestrial realms, which was crucial for our approach.

However, using various sources for the threshold values was inevitable due to the high number of chemicals, as well as different ways authors reported n class="Disease">stressor inpan>tenpan>sity inpan> inpan>dividual studies. We relied on authors’ statemenpan>ts for stanpan>dardization inpan> two mainpan> situations: terrestrial studies that reported pesticide inpan>tenpan>sity inpan> terms of anpan> application rate (n = 6 studies), anpan>d aquatic studies that tested the effect of mixture of a large number of pesticides, anpan>d reported pesticide levels as sums of toxic unpan>its (a method that is similar as what we used to stanpan>dardize pan> class="Disease">stressor levels, where the concentration of each pesticide is divided by standard ecotoxicological data before being summed across pesticides into an overall sum of toxic units) (n = 4). Because application rates and pesticide mixtures are not regulated by our two main authoritative sources, we relied on the authors’ statements (recommended application rates) and expert knowledge (safe levels for toxic units) to standardize stress intensity in those cases.

We agree with the reviewer that the n class="Gene">nutrienpan>t-diversity dataset has a more homogenpan>eous basis for stanpan>dardization compared to the pan> class="Disease">stressor-diversity dataset. However, we would like to stress the fact that our final datasets are all dominated by similar sources for standardizing stressors: thresholds from ECHA or USEPA for 80 and 90% of observations in the stressor-diversity and stressor-abundance datasets, respectively, and for nutrients, thresholds from WFD for 100 and 75% of observations in the nutrient-diversity and nutrient-abundance datasets, respectively. Therefore, it is possible that the results are little affected by the small number of observations for which quality criteria come from a different source than the main ones.

We modified the text to clarify our standardization approach in the Materials and methods section:

“We used consistent sources for the environmental quality criteria as much as possible. […] Despite this, the final datasets were all dominated by similar sources for standardizing stressor anpan>d pan> class="Gene">nutrient intensity levels: thresholds from ECHA or USEPA for 80 and 90% of observations in the stressor-diversity and stressor-abundance datasets, respectively, and for nutrients, thresholds from WFD for 100 and 75% of observations in the nutrient-diversity and nutrient-abundance datasets, respectively.”

We further modified the Discussion section that now discuss the importance of standardization in our results:

“Such contrasting patterns between stressor anpan>d pan> class="Gene">nutrient intensity effects may reflect the greater number of stressor types (different metals, pesticides, mixtures) covered by individual studies compared to nutrients. […] However, future studies focusing on stressor intensity effects across ecosystems would greatly benefit from coordinated efforts to derive quality criteria encompassing the vast and rapidly increasing number of chemical stressors (Wang et al., 2020).”

2) The stressor definpan>ition should be reconsidered. The authors definpan>e: "pan> class="Disease">stressors (e.g., temperature, drought, chemicals) and resource enrichment (e.g., of CO2 or mineral nutrients)." From a theoretical perspective, the temperature would also be a niche-defining resource. From a practical perspective, it is confusing for the reader to contrast "chemical stressors" to CO2 and NO3, the latter of which are clearly chemicals and often anthropogenic. Then both are "environmental change drivers", which makes the terminology even more complicated for the reader. Why not just use stressors and environmental change drivers, which have different niche-defining properties? Some are unimodal (nutrients, temperature), others exhibit a threshold relationship and others may be log-linear.

We agree that the previous wording may have been confusing. Chemical stressors are enpan>vironmenpan>tal chanpan>ge drivers (see e.g. Bernpan>hardt et al., 2017 Frontiers inpan> Ecology anpan>d the Enpan>vironmenpan>t), but contrarily to pan> class="Gene">nutrients, they cannot be consumed. This distinction (being consumed by the focal populations or not) is key to our conceptual framework. Whether or not these variables have implications for niche partitioning is not part of our framework. In that sense, temperature cannot be consumed and therefore does categorize as a stressor in our framework. Therefore, we would like to maintain the adopted terminology. However, we have now changed the wording in the Introduction to define our two main categories in a clearer way, but to also mention potential caveats with any grouping approach. The text now reads:

“We postulate that there are two main categories of environmental change: stressors anpan>d resource shifts. While pan> class="Disease">stressors cannot be consumed, and act as conditions that alter growth rates (e.g., temperature, drought, chemical stressors), resources are by definition consumed (e.g., CO2 or mineral nutrients), which has important implications for how they should enter theory (De Laender, 2018; Chase and Leibold, 2003).”

We further added a statement in the Discussion highlighting that additional groupings of drivers could be addressed by future work focusing on a greater number of drivers:

“Finally, although our comparison of n class="Disease">stressors versus resources allowed us to test a clear concept, anpan>y kinpan>d of groupinpan>g inpan> ecological studies may mask some of the variation withinpan> the categories, anpan>d future studies may be inpan>terested inpan> different categories. […] This would also allow to test additional groupinpan>gs of drivers anpan>d ecological concepts unifyinpan>g pan> class="Disease">stressors and resources (De Laender, 2018; Harley et al., 2017).”