Increasing the number of stressors reduces soil ecosystem services worldwide

Increasing the number of environmental stressors could decrease ecosystem functioning in soils. Yet this relationship has never been globally assessed outside laboratory experiments. Here, using two independent global standardized field surveys, and a range of natural and human factors, we test the relationship between the number of environmental stressors exceeding different critical thresholds and the maintenance of multiple ecosystem services across biomes. Our analysis shows that, multiple stressors, from medium levels (>50%), negatively and significantly correlates with impacts on ecosystem services, and that multiple stressors crossing a high-level critical threshold (over 75% of maximum observed levels), reduces soil biodiversity and functioning globally. The number of environmental stressors >75% threshold was consistently seen as an important predictor of multiple ecosystem services, therefore improving prediction of ecosystem functioning. Our findings highlight the need to reduce the dimensionality of the human footprint on ecosystems to conserve biodiversity and function.

Increasing the number of environmental stressors could decrease ecosystem functioning in soils. Yet this relationship has not been globally assessed outside laboratory experiments. Here, using two independent global standardized field surveys, and a range of natural and human factors, we test the relationship between the number of environmental stressors exceeding different critical thresholds and the maintenance of multiple ecosystem services across biomes. Our analysis shows that having multiple stressors, from medium levels (>50%), negatively and significantly correlates with impacts on ecosystem services and that having multiple stressors crossing a high-level critical threshold (over 75% of maximum observed levels) reduces soil biodiversity and functioning globally. The number of environmental stressors exceeding the >75% threshold was consistently seen as an important predictor of multiple ecosystem services, therefore improving prediction of ecosystem functioning. Our findings highlight the need to reduce the dimensionality of the human footprint on ecosystems to conserve biodiversity and function.
Human activities are leaving a complex footprint of impacts on terrestrial ecosystems associated with multiple environmental stressors (for example, drought, soil acidity or heavy metals), jointly affecting the majority of Earth's ecosystems 1-4 ; yet we know very little about such concurrent effects. This is unfortunate because ecosystems are continuously subjected to a wide range of stressors, including anthropogenic disturbance such as pesticides, microplastics or land-use change 1-7 . A recent systematic mapping revealed that over 98% of papers published on global change and soils, reporting on results from experimental work, covered only one or two global-change stressors 8 . This means that, for a number of reasons, including logistics 9 , the combinatorial explosion problem 10 and an insufficient collaboration among the subfields of global-change biology 1 , experimental work has covered very insufficiently what happens in reality when multiple environmental stressors co-occur and interact. Experimental work covering the effects of multiple stressors on biodiversity-ecosystem functioning has shown that just the number of factors has a large effect 8 . Thus, while single stressors by themselves had relatively minor effects (nominally positive, negative or neutral effects), together negative effects on soil functions were much stronger, that is, synergistic. Working on the Article https://doi.org/10.1038/s41558-023-01627-2 multiple ecosystem functions. Soil ecosystem services are defined from a human perspective and represent a wide range of benefits (carbon storage, food provision and so on). If the number of stressors explains a large portion of variability in biodiversity and process data, this can inform efforts to reduce the dimensionality of human impacts and aid conservation.
In this Article, we explore the relationship between the number of environmental stressors, including a combination of climatic, physico-chemical and anthropogenic stressors, exceeding critical-level thresholds and the maintenance of multiple ecosystem services associated with soil biodiversity and function at a global scale. To provide independent and robust evidence for the link between multiple environmental stressors and ecosystem multiservices, we used two standardized global surveys including a total of 218 ecosystems from six continents and across a wide range of climatic and vegetation types ( Supplementary Fig. 1, Global surveys #1a and #2 and Methods). Our study includes seven groups of climatic (aridity, maximum temperature and seasonality), physico-chemical (salinity, heavy metals and distance from neutral pH) and anthropogenic (human-influence index and fertilization) factors that potentially cause stress (stressors) when passing high levels (determined using a threshold approach as explained in the following). Our Global dataset (Global #1b; n = 48) includes a subset of locations from Global #1a including these seven environmental stressors plus two additional stressors, pesticides and microplastics. We then determined the number of groups of environmental stressors simultaneously passing all possible threshold values (5-95%, based on the maximum observed level for each stressor) following the original threshold approach of ref. 20 . Subsequently, it was possible to test whether ecosystem multiservices are reduced when a number of multiple stressors are above certain values set as thresholds. To assess ecosystem multiservices, we focused on 15 surrogates of ecosystem functions associated with six fundamental ecosystem services (Supplementary Table 2 and Supplementary Figs. 2 and 3): soil biodiversity, organic-matter decomposition, pathogen control, plant productivity, water regulation and nutrient cycling. We characterize multiple ecosystem services using the weighted multiservices 4 approach considering individual organism level, other researchers studying the model plant Arabidopsis obtained similar results 11,12 . This may be the case because with an increasing number of stressors, an increasing proportion of species or component processes may be negatively affected, leading to an overall decrease in both biodiversity and ecosystem process rates. The hypothesis 8 emerging from such experiments is that the number of global-change manifestations, in addition to their actual identity and magnitude, can impact the functioning of terrestrial ecosystems.
Despite the huge importance of understanding how multiple environmental stressors influence the biodiversity and function of terrestrial ecosystems under natural conditions, the relationship between the number of multiple stressors and functioning has not been assessed globally in real-world terrestrial ecosystems across wide gradients of environmental conditions. Testing this hypothesis about stressor number using observational approaches may yield important knowledge on what is the overall actual importance of such stressors in the real world or whether this importance can be masked by natural variability of ecosystems. This knowledge is also critical to understanding the future of terrestrial ecosystems that are suffering from an increasing number of environmental stressors. In addition to experimental approaches, observational field studies can provide further ecological insights as to how multiple stressors influence ecosystem performance. This brings the advantages of increased realism and high external validity, advantages that come at the cost of mechanistic resolution, degree of control and, thus, the ability to demonstrate causality. Observational studies routinely use a wide range of predictor variables to increase the proportion of variance explained or to account for factors that are not of immediate interest but still vary across samples. Studies on soil biodiversity and properties of this nature have in recent years been increasingly carried out at a continental 13,14 or global scale [15][16][17][18][19] , achieving broader insights by covering wide ranges in environmental conditions across biomes. However, such studies have not yet tested, at any scale of observation, the relationship between the number of environmental stressors exceeding critical-level thresholds, including anthropogenically caused influences, and the maintenance of multiple ecosystem services associated with soil biodiversity and confidence intervals of the average standardized coefficients. Global #1a (n = 131 sites) and Global #2 (n = 87 sites) include seven environmental stressors: aridity, temperature, seasonality, salinity, distance from neutral pH, heavy metals and human influence. Global #1b (n = 48) includes a subset of locations from Global #1a with these seven environmental stressors plus pesticides and microplastics.
Article https://doi.org/10.1038/s41558-023-01627-2 six ecosystem services (multiservices). We hypothesized that the number of environmental stressors exceeding medium-high critical threshold levels will influence multiple ecosystem services associated with soil biodiversity and multiple ecosystem functions.

Results and discussion
We first systematically assessed the relationship between the number of stressors exceeding all possible threshold values (5-95%) and ecosystem multiservices (see also Supplementary Figs. 4 and 5 for additional results). Using this approach, we do not presume an a priori level of stress-reducing function (for example, 80%), but calculate the effect of number of stressors on multiple ecosystem services assuming different scenarios iteratively; that is, each threshold is an assumption of the minimum level of stressor required to be detrimental for ecosystem functioning. Our analyses provide new evidence that the number of environmental stressors crossing a high-level threshold ( Fig. 1; see Supplementary Table 3 for P values) is negatively correlated with the capacity of soils to support multiple ecosystem services. Using two independent global surveys, we show that the higher the number of environmental drivers crossing a >50% stressor threshold ( Fig. 1 Table 3 for P values). Our ecosystem services are comparable to those in ref. 8 , suggesting that laboratory experimental outcomes can also be found in terrestrial ecosystems at a global scale under natural conditions. Our framework further expands that in ref. 8 by explicitly considering stressor intensity levels.
Because the effect of multiple stressors could be masked by the intensity of individual stressors, we then used multimodel inference to evaluate whether the number of multiple stressors exceeding different critical thresholds (>25%, >50% and >75%) was always selected among the best models (those with a combination of predictors able to maximize variance absorbed with the minimum set of predictors possible).

Fig. 2 | The number of stressors operating above high levels is an important independent factor for multiservices when compared with individual stressors.
The relative importance of the multistressors index (number of stressors surpassing a given threshold) compared with individual factors to assess multiservices. a, A qualitative assessment of variable importance by showing results from a multimodel inference approach. Lists of the best generalized linear models out of all potential combinations of environmental stressors are given for Global #1a and Global #2. The R 2 s of the saturated models (those from which variable importance was calculated) were 0.3834 for Global #1 and 0.5905 for Global #2. If the variable (columns) was included in the model (is significant), this is shown with a coloured box (green for multistressors, red for climatic individual stressors, blue for soil-related stressors and yellow for other forms of human influence). Effect sizes of each variable can be seen in Supplementary Table 6. The Bayesian information criterion (BIC) indicates the suitability of the model (the lower the better) and delta indicates the difference of BIC with respect to the best model (<4 indicates similar performance to the best model). The weight of the models is also represented in the table and indicates how each of the best models would contribute to an average model (results in Supplementary Table 7). b, A quantitative assessment of variable importance using a variance partitioning scheme. The pie charts indicate the percentage of explained variance of multiservices (R 2 ) that is associated with each predictor.
Article https://doi.org/10.1038/s41558-023-01627-2 This analysis tested more thoroughly whether the importance of the number of multiple stressors exceeding a threshold level was consistent even after controlling for individual stressors. We conducted these analyses in Global surveys #1a and #2, supporting the largest number of sites. We found that the number of stressors exceeding the >75% threshold together with human influence were the most consistently selected factors in all the best models predicting the distribution of multiple ecosystem services in two independent global surveys ( Fig. 2). Heavy metals were also selected in many of the best models from both global surveys, highlighting their importance for explaining multiservices (Fig. 2). The number of stressors exceeding the >75% threshold (also >25 and >50%, depending on the dataset examined) was maintained in the models even when enforcing individual stressors to be kept in the model (Supplementary Table 7). In general, multiple individual ecosystem services such as the biodiversity of soil tardigrades, Platyhelminthes and mycorrhizal fungi, potential infiltration rates, net primary productivity and soil-borne pathogen control had consistent responses to multiple environmental stressors, declining in response to a high number of stressors; however, other individual services gave rise to database-dependent responses ( Fig. 3 and Supplementary Fig. 3). Thus, increases in the number of stressors had different influence on the diversity of annelids, Arthropoda and fungal decomposers and glucose and lignin respiration in Global Survey #1 versus Global Survey #2 (Fig. 3). The negative influence of multiple stressors on the diversity of nematodes was especially important when considering stressors such as microplastics and pesticides in Global Survey #1b. This may be related to the fact that Global Survey #1 includes multiple urban areas, which are known to be exposed to multiple additional environmental stresses, many of which have not been measured here; thus, some of these unmeasured stressors may have contributed to this observed effect. Also of interest was the consistent negative response of potential infiltration to multiple stressors in Global Surveys #1 and #2, even when these stressors occurred at a relatively low level of intensity (>25%; Fig. 3). Finally, the number of multiple stressors was consistently positively correlated with the availability of phosphorus in soil, which is probably associated with human-driven fertilization processes.
With the analyses performed so far, we could investigate the overall importance for ecosystem services of the number of stressors crossing a given threshold. However, we also wanted to explore the possibility that the effects of multiple stressors could depend on the individual climatic stressors. Thus, the response of multiple ecosystem services to climate change (individual climatic stressors) could be additionally impacted (if effect is additive) by multiple environmental stressors, including human multidimensional disturbances, or directly modulated by them (if effect is interactive). We built linear models using as the response variable multiple ecosystem services and, as predictors, climatic groups of variables and their interaction with the multistressor index at the 75% threshold (MS75). We chose MS75 because this threshold showed a clear effect in the multimodel inference approach compared with MS25 and MS50 (Fig. 2). Apart from assessing the significance of such interactive effects (Supplementary Table 6), we also plotted them for visual inspection. To do so, we divided the data according to the number of functions surpassing 75% of their maximum and provided a plot of the ones with a high number of functions surpassing this threshold and another of those with a low number of functions surpassing this threshold (Fig. 4). We found that, in general, in both global surveys, increases in maximum temperature and aridity can have important additive negative influences on the maintenance of multiple ecosystem services in terrestrial ecosystems when interacting with multiple environmental stressors, including human influence and soil stressors. This pattern was not observed, however, for maximum temperature in Global Survey #2. This finding supports the notion that climate change impacts on soil biodiversity and function could be larger than previously expected as there are additional pressures from other human  influences and abiotic factors that need to be simultaneously considered when evaluating the impacts of climate change on biodiversity and ecosystem services. Our global study presents an important advance in current knowledge on the role of multistressors in driving soil ecosystem services by explicitly considering the concepts of number of stressors and stressor intensity in real-world ecosystems. Previous experimental work inevitably included a limited suite of factors (up to ten for soil microcosms 8 and up to six for a model plant system 11 , with two fixed levels of each factor). More complexity is difficult to achieve within the confines of experimental design and logistics, and the necessary number of replicates, at least in ecology. By contrast, the work presented here fundamentally reflects a naturally occurring continuity of factor levels that we converted into numbers of factors by means of a threshold approach. In addition, all other influences were in effect, and the study thus spans a very wide range of climate and soil types, as well as vegetation feedbacks. In addition, by using two independent global surveys, we provide robust evidence that the number of multiple stressors can negatively influence multiple ecosystem services across environmental gradients. Finding a similar pattern for both surveys lends substantial additional credibility to the 'number of factors' hypothesis.
It will be important to develop policy instruments or approaches that take into account the overall dimensionality of human-caused factors 8 . Currently, much policy work focuses on individual factor types, for example, microplastics (for example, EU-level initiatives to restrict microplastic pollution) or the complex suite of climate impacts including drought and warming (Intergovernmental Panel on Climate Change, https://www.ipcc.ch/). Our findings indicate that a growing number of human-induced stress factors such as pesticides, microplastics and general human influence/disturbance can further exacerbate the negative impacts of climate change on ecosystem performance. A recent study demonstrated that fungicide application impairs the ability of beneficial soil fungi to supply nutrients to their host plants 21 . It is important to assess whether such effects would be even stronger when these fungi are simultaneously exposed to drought, microplastics, heavy metals and other soil stressors. Moreover, our findings indicate that ecosystems supporting a high number of environmental stressors passing high levels of stress could result in important declines in biodiversity and function. Thus, our findings strongly suggest that to protect and conserve terrestrial ecosystems and their contributions to people, there needs to be a shift in focus towards entire suites of influences, in terms of both work conducted in environmental sciences and relevant policy regulations. Our study suggests that the number of stressors largely shapes the soil multiservices worldwide, which needs to be considered in future soil macroecology studies. This knowledge is critical to anticipate reductions in soil biodiversity and function in response to global change and across terrestrial ecosystems, which could be overlooked when global change drivers are considered individually.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-023-01627-2.

Methods
The analyses in this study 22 are based on two independent global field surveys: Global Survey #1 and Global Survey #2. These databases were selected because they include a wide range of similar soil biodiversity and ecosystem functions measured in the same samples. This is critical when exploring multifunctionality indices and not very common in the global-scale literature.

Global Survey #1
Three composite soil samples from five soil cores (top 5 cm) were collected from the most common microsites (under vegetation and in open areas between plant canopies) in 131 ecosystems across 6 continents ( Supplementary Fig. 1) ranging from deserts and polar areas to tropical and temperate forests. This sampling was conducted between 2017 and 2019. This survey targeted 59 well-established urban parks and large residential gardens (urban forests and lawns) and 72 relatively undisturbed natural ecosystems. A subset of the locations in this survey has been previously used to investigate global homogenization patterns in soil biodiversity 23 .

Global Survey #2
Five composite soil samples from five soil cores (top 10 cm) were collected from the most common microsites (under vegetation and in open areas between plant canopies) in 87 terrestrial ecosystems located in 9 countries on 6 continents ranging from deserts and polar areas to tropical and temperate forests ( Supplementary Fig. 1). This sampling was conducted between 2016 and 2017. A subset of the locations in this global survey is available from the literature and has been previously used to investigate linkages between soil biodiversity and ecosystem function 4 .

Ecosystem services and functions
Our study includes 15 ecosystem variables included in 6 ecosystem services: organic-matter decomposition (soil respiration and lignin-and glucose-induced respiration), soil biodiversity (diversity of annelids, arthropods, nematodes, Platyhelminthes, tardigrades, mycorrhizal fungi and fungal decomposers), pathogen control, plant productivity (net primary productivity), water regulation (potential infiltration) and nutrient cycling (available nitrogen (N) and phosphorus (P)) (Supplementary Table 2). These ecosystem services were selected for two reasons: (1) they constitute good surrogates of multiple ecosystem services associated with plant production, organic-matter decomposition and soil biodiversity, and (2) they are directly comparable to services in ref. 8 , for which experimental evidence exists on the role of multiple stressors in driving soil biodiversity and function. Soil respiration rates (basal respiration and glucose-and lignin-induced respiration) were determined using a composite sample per plot using Microresp at 25 °C and 60% water-holding capacity. The biodiversity of soil invertebrates (richness of soil invertebrates; amplicon sequence variants (ASVs)) was obtained from all composite soil samples using Miseq Illumina sequencing analyses by characterizing a portion of the eukaryotic (invertebrates) 18S ribosomal RNA genes using the Euk1391f/EukBr (https://earthmicrobiome.org/ protocols-and-standards/18s/) primer sets. The diversity of mycorrhizal fungi (arbuscular and ectomycorrhizal fungi) and fungal decomposers was obtained from PacBio sequencing analyses (all soil samples for Global Survey #1 and a composite sample per plot for Global Survey #2). Molecular analysis of the full-length internal transcribed spacer (ITS) region for fungi was performed using ITS9mun/ ITS4ngsUni primer 24 . Bioinformatic processing was performed using USEARCH 25 and UNOISE3 26 . Phylotypes (ASVs) were identified at the 100% identity level. The ASV abundance tables were rarefied to 5,000 (18S) and 1,000 (ITS) sequences per sample. FungalTraits was used to characterize different fungal functional groups (decomposers, mycorrhizal and potential plant pathogens) in our ASV relative abundance ITS table. Pathogen control (inverse of proportion of fungal plant pathogen as defined in ref. 4 ) was also determined from the PacBio ITS ASV table 27 . PacBio technology supports enough resolution to determine potential soil-borne plant pathogens. Net primary productivity was estimated from the normalized difference vegetation index using satellite imagery (Moderate-Resolution Imaging Spectroradiometer) at 250 m resolution. Potential water infiltration was determined in the lab in a composite sample per plot using a method similar to that described in ref. 4 . In Global Survey #1, available N and P were estimated in all composite soil samples using root simulators in the lab and a mix of soil and water (1 h extraction). In the case of Global Survey #2, a salt-extraction (potassium sulfate 0.5 M) was used to extract soil inorganic N, and sodium bicarbonate was used to extract inorganic P. In both cases, soil N and P were determined colourimetrically using the methods described in ref. 4 .

Plot estimates of environmental factors and ecosystem services
Before statistical analyses, within-plot information on all environmental factors (for example, heavy metals, pH and salinity), functions and soil biodiversity (richness of soil organisms) metrics, derived from three (Global Survey #1) or five (Global Survey #2) composite soil samples per plot, was averaged to obtain plot-level estimates. This approach allowed us to work at the site level, matching with environmental factors such as climate and human influence. This was not needed for those analyses including a single composite soil sample per plot.

Weighted ecosystem multiservices
We used average weighted ecosystem multiservices. To determine this measurement of multiservices, we first standardized our functions between 0 and 1 as done in ref. 28 . We averaged all functions within six ecosystem services and then calculated the average of six services as our metric of multiservices. By doing so, we aim to ensure that all services equally contribute to our multiservices index and that the number of functions within each service is not influencing our results.

Groups of natural and human-based environmental stressors
In this study, we worked with seven groups of stressors (instead of multiple individual factors within each group) for two reasons. First, individual factors within each group (for example, within heavy metals) are highly correlated with each other, suffering from multicollinearity. Second, different groups of environmental stressors comprise a different number of individual factors (for example, from one to six). These seven groups of stressors were selected on the basis of two criteria: their well-known importance and data availability in the used databases. Our groups of seven stressors do not suffer from statistical multicollinearity presenting different types of stressors, and they reflect largely independent statistical entities (Supplementary Table 4 and Supplementary Figs. 2 and 3). In addition, working at the individual factor level would not allow a fair and weighted assessment of the influence of the number of environmental stressors exceeding different stressor thresholds, which could result in groups having more individual variables also giving rise to larger and overrepresented impacts. The selected groups of stressors showed relatively low levels of correlation, suggesting that they represent different dimensions of environmental impact (Supplementary Fig. 2).
We considered seven common environmental factors that can result in environmental stress when passing high levels: aridity (inverse of aridity index), temperature (maximum temperature), seasonality (precipitation and temperature seasonality and diurnal temperature range), salinity, distance from neutral pH, levels of heavy metals (soil As, Cd, Cr, Cu, Ni, Pb and Zn) and human influence (for example, human-influence index and fertilization). See Supplementary Table 1 for details. These data were available for Global Survey #1a and Global Article https://doi.org/10.1038/s41558-023-01627-2 Survey #2 (Supplementary Fig. 1). We also considered level of pesticides and microplastics, for which data were available for a subset of locations in Global Survey #1 (Global Survey #1b; Supplementary  Fig. 1 and Supplementary Table 1). This information was included in our manuscript to account for the potential influence of well-known stressors such as pesticides and microplastics when data were available.
In both global surveys, aridity Index was obtained from ref. 29 (v.2). Temperature and seasonality variables were obtained from WorldClim (v.2) (http://www.worldclim.com/version2). Electrical conductivity was measured in all soil samples as a surrogate of soil salinity as done in ref. 30 . Soil pH was measured with a pH meter in a soil/water solution (1/10). The concentration of As, Cd, Cr, Cu, Ni, Pb and Zn was measured as surrogates of heavy metals using inductively coupled plasma optical-emission spectrometry with Thermo ICP 6500 Duo equipment (Thermo Fisher Scientific) before a microwave digestion. In the case of Global Survey #1, human influence includes information on city influence (natural = 0; urban = 1), fertilization (0 = no; 1 = yes) and mowing (0 = no; 1 = yes) collected in the field. In the case of Global Survey #2, we collected information on nitrogen fertilization in manure 31 and via direct human application 32 . For both databases, we also considered the human-influence index 33 . This index aggregates important information on eight factors: major roadways, navigable waterways, railways, crop lands, pasture lands, the built environment, light pollution and human population density. The concentration of microplastics was measured using RAMAN spectroscopy ( Jasco NRS-5100) with laser excitation frequency and intensity 784.79 nm and 11.8 mW, respectively, following a density extraction of the microparticles with a saturated NaCl solution. The RAMAN spectra of the analysed particles were compared with reference polymers from the spectral library Open Specy 34 .

Number of environmental stressors over a threshold
We then calculated the number of stressors passing three thresholds (for example, 25%, 50% and 75%) of maximum stressor levels (see ref. 21 for a similar approach in a different context). This approach assumes that when a high number of environmental groups of factors (for example, heavy metals and seasonality) pass high levels (for example, >75%; versus maximum values), this can have negative impacts on both biodiversity and function. Before we determined the number of stressors passing over these three thresholds, we standardized all individual environmental stressors between 0 and 1 (within each dataset) and averaged all environmental stressors (for example, maximum temperature) within each group of stressors (that is, temperature) (Supplementary Table 1). By doing so, we aimed to have different groups of stressors equally contributing to our multistressor indexes, avoiding an over-influence of stressor groups with multiple environmental variables. We then determined the number of stressors exceeding the thresholds for each location (5-95%), as done in ref. 20 , and used this number for further analyses as the explanatory variable. For Fig. 1, we used the same approach, but we did not restrict the number to 25%, 50% and 75% thresholds; we spanned the range of 1%-99% thresholds in unitary steps. Our analyses are based on seven measured environmental stressors in the cases of Global #1a and Global #2 and nine environmental stressors measured in the case of Global #1b (Supplementary Table 1).

Statistical analyses
We first used linear regressions and Spearman correlations to evaluate the link between the number of stressors passing the three thresholds and ecosystem multiservices, the 6 ecosystem services and the 15 ecosystem functions.
To assess the relationship of multiple stressors acting together on ecosystem multiservices, we correlated the number of stressors passing a given threshold value (from 5% to 95% iteratively increasing by 1%) with the multiservices index. This yielded a total of 91 regressions, each associated with a threshold value. Each of these regressions was bootstrapped 100 times, and the slopes were obtained to be plotted against the threshold for testing their significance (when confidence intervals at 95% of the bootstrapped slopes do not intersect the value of 0, they are significant; see Fig. 1b).
By increasing the threshold level in this analysis, we were able to see whether multiservices respond mostly to critical levels of stressors. For example, a relationship turning significant when a threshold of 20% is reached means that the variance of number of stressors performing at 20% of their standardized values can exert a significant effect on ecosystem functioning. Moreover, this stress level is not associated with one stressor in particular, but with several acting synergistically (one stressor can have high values and the multistressor index be very low if other stressors are mild).
Because significance of a multistressor indicator does not necessarily imply that the multistressor is a better predictor than the isolated stressors, we then conducted a multiple linear regression using all stressors in isolation and the multistressor index for 25, 50 and 75% thresholds for the database of Global #1a and Global #2. Analysis of variance inflation factor (VIF < 5) revealed a lack of multicollinearity between these predictors, which was also confirmed by overall low correlations between the individual stressors (Supplementary Figs. 2 and 3 and Supplementary Table 5). The (saturated) model for Global #2 had a term with VIF > 5, namely, M25; we thus reran the model, excluding this term, and still found the overall result to hold. To test the importance of the multistressor indices, we performed two evaluations using this multiple regression as a baseline.
First, we measured the proportion of variance explained by each of the predictors (multistressor indices plus individual stressors) using a variance partitioning scheme with the package relaimpo 35 in R. The results of this analysis inform about the corresponding proportion of the variance explained by each predictor of the model. Usually, predictors absorbing large amounts of variance from the response variables are considered more important.
Second, we performed a multimodel inference diagnosis using the package MuMIN 36 in R. This analysis conducts model fitting of all possible combinations of predictors in a given multiple regression model and extracts the Bayesian information criterion (BIC) of each combination to be used for ranking the models according to the principle of maximum parsimony. BIC informs about a trade-off on the variance explained by the model and the number of parameters required for its fitting, penalizing models that are too complex. The best model (lowest BIC) is ranked first, and the rest of the models are ordered on the basis of the deviation of their BIC. In general, models with a difference in BIC lower than 2 are considered as having similar performance. To use a more conservative criterion, we focused on models with BIC difference up to 4. The examination of which variables are included within the selected best models informs about the importance of such variables. This approach is complementary to variance partitioning because it is less dependent on the amount of variance absorbed and more on whether this variance is absorbed exclusively by a particular predictor (an important variable, able to explain a unique proportion of the variance, even if that variance amount is very little, is usually not dropped from the best models).
We additionally performed a canonical correspondence analysis (using vegan 37 R package) that confronted the values of our individual functions with contrasting levels of multistressor thresholds (from 5 to 95 taken each 5%). Canonical correspondence analysis is able to order the six functions measured in our study in a canonical two-dimensional space characterized by the combination of multistressor levels that best separates functions from each other. The result can represent functional trade-offs emerging when considering the linear effects of multistressors on the functions (for example, separating functions that respond mostly to high thresholds of multistressors from functions that respond to lower thresholds).
All analyses were done with R 38 .

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
The database used in this manuscript is available at Figshare.com (https://figshare.com/s/e58852c4ffe5c3553a42) 22 .

Code availability
The R codes used in this manuscript are available at Figshare.com (https://figshare.com/s/e58852c4ffe5c3553a42) 22 .