Literature data on natural plant litter decomposition

We conducted an extensive literature survey for peer-reviewed publications in the ISI Web of Science database (http://apps.webofknowledge.com/) on September 1st 2023. We used (warming OR heat* OR OTC OR open top chamber*) AND (litter* OR litter bag) AND (decomposition OR mass loss) as search criteria, which returned 1184 studies (Figure S1). We considered terrestrial field studies that compared litter decomposition (mass loss and decomposition rate) under  experimentally increased temperatures (methods found in our search were open-top chambers, heating cables, infrared heaters, sunlit controlled-environment chambers, UVB filter films, open-topped polythene tents, and closed-top chambers) and ambient conditions. From the 60 studies that met our criteria, we extracted mean values, sample sizes and measures of variation (i.e., standard errors or standard deviations) for decomposition (i.e., decomposition rates, absolute and relative mass loss, remaining mass of plant material). We contacted the corresponding authors to obtain access to the raw data for studies that did not report them and had to exclude 8 studies due to insufficient reporting. This resulted in 52 studies used for the meta-analysis (Table S1). Whenever warming was applied in factorial combination with one or more additional treatments (e.g., warming and plant species removal), we extracted the parameters of interest for the warming treatment only, together with the ambient control. If the litter was incubated at different time steps, each time step was used as an independent data point. We thus extracted a total of 523 paired (ambient vs warmed) data points from the 52 studies, either directly from the text or tables or from figures using the software WebPlotDigitizer (v. 4.6, Rohatgi 2021). When decomposition was reported as remaining mass of plant material, the latter was transformed into mass loss.

We extracted coordinates of each study location (Figure 1A), the incubation duration of the litter (from 14 days to 4.9 years, standardised to days), the mesh size of the litter bags (from 0.02 to 5 mm), the position of incubation (i.e., if litter bags were put on the soil surface or buried below ground), the plant species and thus the plant functional type (i.e., forb, nonvascular, graminoid, woody species), and the plant organ type (i.e., leaf, shoot or root). For 32 studies, we also extracted litter C:N ratio reported by the researchers (ranging from 12 to 201). All reported values were from single species, with the only exception being two studies on root decomposition, which included a mixture of grass species. Yet, as all the species in these samples were indeed within the graminoid functional type, we included these studies in the meta-analysis. For each study, we extracted the duration of the warming experiment prior to incubation start (from first year to 23 years). The warming method was classified as Heating cables (number of studies n=11), Infrared heaters (n=17), and Open-top chambers (n=19), with ‘Other methods’ including Sunlit controlled-environment chambers (n=1), UVB filter films (n=1), Open-topped polythene tents (n=2) and Closed-top chambers (n=1).

Standardised plant litter data from open-top chamber warming experiments

Following the standard Tea Bag Index protocol (Keuskamp et al., 2013), green (Camellia sinensis; EAN no.: 8 722700 055525) and rooibos (Aspalathus linearis; EAN no.: 8 722700 188438, Lipton, Unilever) tea bags with woven nylon mesh (0.257 mm), were buried at a depth of 8 cm and at a distance of at least 15 cm from each other in open-top chambers (OTC) and controls at 57 locations (Figure 1, Table S2). The incubations covered one growing season (82 ± 18 days; mean ± SD), that is, from May/June 2016 to August/September 2016 in the northern hemisphere and from January 2017 to March 2017 in the southern hemisphere. For two sites in Japan (i.e., JPN_1 and JPN_3, Table S2), tea bags were incubated from July to October 2012. Retrieved bags were cleaned of adhering soil and roots. The mass of the remaining tea was determined after drying it in an oven at 60-70 °C for at least 48 h. To align with the literature data, we calculated treatment means of mass loss, sample sizes and standard deviations for each experiment/GPS location.

Explanatory macro-environmental drivers

We obtained map-based environmental data based on the geographical locations of the study sites to identify macro-environmental factors that may influence the response of decomposition to warming. We used 48 environmental layers reflecting major gradients in climate, soil, vegetation, and topographic variables as covariates in our analysis (Table S3).

Due to the occurrence of many confounding variables, we summarised the macro-environmental variation across studies with a Principal Component Analysis (PCA; Table S3). The first principal component (PC 1) was strongly positively correlated with temperature-associated variables and negatively correlated with soil organic carbon (SOC) and explained 26.9 % of the total variance (Figure 3A, Table S3). The second component (PC 2) correlated positively with precipitation-associated variables and explained 18.1% of the total variance (Figure 3A, Table S3). The third PC axis was not considered as it described negligible amounts of the variation (4.2%). In our dataset, the range of annual mean temperature was -12 to 28 °C, annual precipitation was 78 to 2100 mm, and soil saturated water content was 42 to 81 %. 

We created four ‘macro-environmental classes’ based on the origin of the PC1 and PC2 variables as a separation line. These four ‘macro-environmental classes’ were described as following: (1) high temperatures and high precipitation (number of effect sizes k=156), (2) high temperatures and low precipitation (k=170), (3) low temperatures and high precipitation (k=156), and (4) low temperatures and low precipitation (k=155) (Figure S2, Table S4).

Explanatory micro-environmental drivers altered by experimental warming

For both datasets (i.e., natural and standardised plant litter), we collected available data on the actual degree of warming, i.e., the mean absolute temperature difference between the warmed and ambient control, as well as soil moisture in warming and control treatments when available. The degree of warming included air or soil temperature measures, depending on whether the litter was incubated on the soil surface or below ground, respectively. We calculated relative change in soil moisture with warming according to:
"relative change in soil moisture  = " ("Mc" /"MW" -1)×100
where MC and MW are soil moisture in control and warming treatment, respectively. Positive and negative values indicate drier and wetter conditions under warming than under ambient conditions, respectively.

Litter quality

We focused on three different, frequently used characterisations of litter qualities: the C:N ratio before decomposition, reported in the initial studies, the decomposability measured as decomposition rate under ambient condition (i.e., standardised to mass loss in % d-1) (Cornelissen et al. 2004; Freschet et al. 2012), and plant functional type (Dorrepaal et al. 2005). We categorised the plant species into four different plant functional types (sensu Chapin et al. 1996), forbs (number of studies n=7), graminoids (i.e., grasses and sedges, n=28), woody species (i.e., shrubs and needle-leaved and broad-leaved trees, n=27), and nonvascular (i.e. mosses, n=4; lichens, n=1). For graminoids and woody species, we were able to further specify litter type into aboveground (i.e., shoots and leaves of graminoids, n=25; broadleaves and needles of woody species, n=25) and below ground plant organs (i.e., roots of graminoids, n=6; and root of woody species, n=2).

Data analysis

In order to evaluate the relative effect of experimental warming on decomposition, we used Hedges' g, which is a standardised mean difference (SMD) calculated by dividing the difference between the mean mass loss in the warming treatment (x ̅1) and control (x ̅2) by the pooled standard deviation (Hedges 1981). 
Hedge's g =(〖(x ̅〗_1-x ̅_2))/(√((n_1-1)*s_1^2  + (n_2-1)*s_2^2) / (n_1+n_2-2))
where n1 and n2 are sample size, and s12 and s22 are the sample variance of the warming treatment and the control, respectively. Therefore, a SMD larger than one indicates that warming enhanced decomposition, while a SMD lower than one indicate that warming decreased decomposition. By using the SMD as a measure of effect size, we were able to synthesize data measured on different scales or units (e.g. mass loss vs. decomposition rate), while still accounting for the precision (variance) of the measurement. We used the escalc() function in the R package METAFOR (Viechtbauer 2010) to derive SMDs and corresponding 95% confidence intervals (CI) from all paired (i.e., warmed/control) mean values, sample sizes and standard deviations for decomposition. We used multivariate linear mixed-effects models fitted through the rma.mv() function (METAFOR package) to calculate the pooled average SMD across all studies. In our multivariate linear mixed-effects models, we used a sampling error covariance matrix to account for the correlation between sampling errors within studies from which multiple effect sizes were extracted. To test whether decomposition was significantly affected by warming, we considered a pooled average effect size estimate to be significantly different from zero if the 95% CI around the mean did not include zero. Effects were considered significant at α=0.05. 
To account for spatial autocorrelation between the study locations, we used a random effect consisting of longitude and latitude (in decimal degrees, with negative values for West and South) based on their great-circle distance (WGS84 ellipsoid method). In multivariate linear mixed-effects models, total heterogeneity (QT) in effect sizes can be partitioned into heterogeneity explained by the model structure (QM) and unexplained heterogeneity (QE). We used the Test of Moderators (QM test) to determine whether there were any significant effects of various moderators  on decomposition SMD (Koricheva et al. 2013).

We tested first for differences between the natural and the standardised plant litter dataset by using the data type (i.e., natural litter or standardised plant litter) as a moderator in the multivariate linear mixed-effects model. Because the effect of warming on decomposition of natural and standardised plant litter did not significantly differ (moderators' test: QM (df = 2) = 2.7, p = 0.26), we combined   the natural and standardised plant litter dataset in the following analyses and added the dataset type to the random-effects structure of our models to account for the different origins of our datasets (random = ~ LON + LAT | dataset type).

To test the impact of macro-environment on the warming effect on decomposition, we first used multivariate linear mixed effects models (n=48) to explore whether the macro-environmental factors individually had a significant effect on the decomposition SMD (Table S6). However, as most environmental factors were confounded, we combined the macro-environmental factors to the underlying gradients using a Principal Component Analysis (PCA) on the scaled environmental variables using the R package FACTOMINER (v.2.4; Lê et al. 2008). We then used the four ‘macro-environmental classes’ created based on the origin of the PC1 and PC2 variables as a separation line, as moderator in the following multivariate linear mixed effects models to test whether the four environmental classes differed in their warming effect on decomposition. We used this factor ‘class’ as interacting moderator in the model totest for interactions in the macro-environment and the natural and standardised plant litter dataset. 
We tested the impact of experimental induced changes in micro-environment, i.e., degree of warming and warming-induced changes in soil moisture and their interaction by testing them as moderators in the same multivariate linear mixed effects model (METAFOR package). We included the ‘macro-environmental class’ as interacting moderator to the model to test whether experimental induced changes in micro-environment differ between the four macro-environmental classes. To investigate whether experimental warming affected temperatures and soil moisture, we used a one-sample t-test to test whether the absolute difference between the warming treatment and the ambient control differed significantly from zero.

To test for differences in the warming effect between the different warming methods used in the different studies and experiments (Table S1, 2), we used ‘warming method’ as moderator in another multivariate linear mixed effects model. In this model, the macro-environmental class was not integrated because the warming methods were not evenly distributed across the four macro-environmental classes (e.g., more OTC studies in higher latitudes). To test for differences in the warming methods in their effect on micro-environment, we used linear mixed-effects models (R package LMERTEST, v. 3.1-3; Kuznetsova et al., 2017) to test the overall effect of the categorical independent variable ‘warming method’ on the continuous dependent variables ‘degree of warming’ and ‘warming-induced changes in soil moisture’, respectively. We used Tukey HSD post-hoc tests (R packages MULTCOMP, v. 1.4-19; Hothorn et al. 2008, and EMMEANS, v. 1.7.5; Lenth 2019) to check for significant differences between the warming methods in degree of warming and warming-induced changes in soil moisture, respectively. We further tested with a linear regression for correlations between warming-induced changes in soil moisture and the degree of warming.

To test for differences in litter quality, described as C:N ratio or ambient decomposability, between plant functional types in different macro-environments, we used linear mixed-effects models (R package LMERTEST) to test the overall effect of the categorical independent variables ‘plant functional type’ (including plant organ types) and ‘macro-environmental class’ and their interactions on the continuous dependent variables ‘C:N ratio’ or ‘ambient decomposability’ (different model for each). We then used Tukey HSD post-hoc tests (R packages MULTCOMP, v. 1.4-19; and EMMEANS, v. 1.7.5) to check for significant differences between the plant functional types and four macro-environmental classes in C:N ratio and ambient decomposability, respectively.

To test our hypothesis that lower litter quality is associated with a stronger positive warming effect on decomposition, we used multivariate linear mixed-effects models (METAFOR package) with the three applied proxies for litter quality ‘C:N ratio’, ‘ambient decomposability’ and ‘plant functional type’ (different model for each of the proxies) as moderators and ‘macro-environmental class’ again as an interactive factor.

In addition, we tested the site-specific drivers related to environmental conditions (absolute latitude and, altitude), experimental setup (duration of warming before the experiment, mesh size) as individual moderators fitting separate multivariate linear mixed-effects models (Table S5). 
For each model, we tested the assumptions of normality and homogeneity of variance of the residuals. Whenever necessary, data were log-transformed (C:N ratio) or rank-transformed (warming-induced changes in soil moisture), with the latter case resulting in a non-parametric regression. Graphical displays were produced using the R packages GGPLOT2 (v. 3.3.6, Wickham (2016)) and ORCHARD (v.2.0, Nakagawa et al. 2021). We used R version 4.2.3 (R Core Team 2023) for all analyses.