We primarily relied on the “tidyverse”61, “metafor”30, “lmerTest”62, “effects”63, “maptools”64, “sf”65, and “raster”66 packages of the R software67 for organizing and analyzing data, and visualizing the results.
To support the use of meta-analysis to quantify the effects of species richness on plant litter decomposition (i.e., species-mixing effects), we searched the literature using the ISI Web of Science (WoS) database (https://clarivate.com/products/web-of-science/; up to and including December 2018). We used a combination of “decomposition” AND “litter”. These keywords matched 12,278 publications. We then reduced this list of the literature with the keywords of “mix* litter,” OR “litter diversity,” OR “litter mix*,” OR “litter species diversity,” OR “litter species richness,” OR “species mix*,” OR “mix* plant litter,” OR “multi* species litter,” OR “mixing effect*,” resulting in 416 publications. We also searched the literature using Scopus (https://www.scopus.com/home.uri) and Google Scholar (https://scholar.google.ca/) and the latter set of keywords, resulting in 765 publications. An anonymous expert identified 15 additional publications. We decided not to include reports in the grey literature among the papers we found, as we could not verify whether, as in all WoS-listed papers, the papers had undergone independent peer review as a quality check. We read through the papers carefully to select those that focused on quantifying the effects of litter diversity on decomposition rates based on a litter-bag experiment using both mixed- and mono-species litter. We focused on the publications that reported values of either mass loss (or mass remaining) or a decomposition rate constant k that quantified the decomposition rate. Specifically, some papers for litter-mixing effects reported only for the additive effects of mixture and had no reports on mass loss or the constant k. In some cases, no sample size nor standard deviations/errors were reported. Mass loss or standard deviations/errors values were reported only for mono-species or multi-species litter. These publications were not included in the present analyses. In cases that had no report on means but instead showed median values, we estimated means and standard errors based on individual data points or percentile values of boxplots. Note that, standard deviations/errors were often invisible because they were too small to read and thus hidden by their data points; in this case, a conservative approach was adopted by using possible maximum values of deviations/errors (i.e., the size of these symbols). Based on these criteria, we extracted information required for our meta-analysis (see below) from a total of 151 studies (Supplementary Data 1). Out of these 151 studies, 131 and 45 studies measured mass loss and the decomposition rate constant k, respectively. Also see the PRISMA work flow diagram (Supplementary Fig. 6).
For the selected studies, we identified the sets of comparisons between mixed- and mono-species litter bags in each publication. That is, a single study could have multiple treatments that focused on different litter substrates, using different mesh sizes for litter bags, having multiple experiments in different locations, changing abiotic conditions, and so on. The decomposition rate should be comparable in a given treatment under the same set of environmental conditions except for differences caused by a different number of plant species. Even for a given treatment, decomposition rates were often reported multiple times in the literature; in such cases, we only included comparisons for litter retrieved at the same time (i.e., after the same incubation period). Note that in a given treatment, the same mono-species litter could be used for multiple comparisons; for instance, mass loss from a litter bag that contained a three-species mixture (species A, B, and C simultaneously) can be compared with the mass loss in three different mono-species litter bags (only species A, B, or C). This could lead to issues related to pseudo replication, which we considered based on a multilevel random effects meta-analysis (see Data analyses). As a result, we identified 6535 comparisons (across 1949 treatments from the 131 studies) for the values of mass loss and 1423 comparisons (across 504 treatments from the 45 studies) for values of the decomposition rate constant k.
After identifying the sets of comparisons, we extracted the sample size (number of litter-bag replicates, n), and the mean and standard deviation (SD) of the decomposition rate from the main text, from any tables and figures, and from the supplemental materials of the selected studies. If standard errors or 95% confidence intervals (CI) were given, we transformed them into SD values. If only figures were given, we used version 3.4 of the Webplot-digitizer software (https://automeris.io/WebPlotDigitizer/) to extract these parameters from the graphs. We also recorded the study location (longitude, latitude, elevation), climatic region (subarctic, boreal, temperate, subtropical, tropical, or other), ecosystem type (forest, grassland, shrubland, desert, wetland, stream, seagrass, lake, or ex situ microcosm), and litter substrate (leaf, root, stem, branch, straw, or other). Microcosm studies were further divided into two categories: terrestrial or aquatic. The former and latter, respectively, placed litter bags on the soil and in water. For litter substrate, most studies used terrestrial plants, but two studies used macrophytes for their decomposition experiment; these data were included in the main analysis and excluded from the subsequent analyses that focused on leaf litter. Note that authors reported the decomposition rate for mixed-species litter bags in different ways, as the values for all species together or for individual species in the same bag. Removing data based on the latter classification had little effect on our results, so we retained that data. Because of the limited data availability, we did not consider the potential influences of species richness (here, the number of litter species; more than 74% of the comparisons used two- or three-species mixtures) and instead considered species richness as a random term (see Data analyses). Based on the geographic locations of the studies, we estimated their present bioclimatic conditions based on the WorldClim database (www.worldclim.org).
Species-mixing effects on litter decomposition
We calculated the unbiased standardized mean difference (Hedges’ d)30 of the decomposition rate between the mean values for the mixed- and mono-species litter. Hedges’ d is a bias-corrected, unit-free index that expresses the magnitude of the deviation from no difference in the response variable between comparisons. Note that many studies performed multiple comparisons for the decomposition rate between mixed- and mono-species litter, potentially causing issues of pseudo-replication and non-independence in the dataset, as is often the case for ecological syntheses68. We thus applied a multilevel random effects meta-analysis30,69 to account for this problem. In a random effects model, the effect sizes for individual comparisons are weighted by the sum of two values: the inverse of the within-study variance and the between-study variance. The multilevel model can also account for a nested structure in the dataset (in which different treatments and multiple comparisons were nested within each study) and is thus appropriate for dealing with non-independence within a dataset. To calculate the effect sizes, we defined their values to be positive for comparisons in which decomposition was faster in mixed-species litter than in mono-species litter and negative when mono-species litter decomposed faster. This is based on a species gain perspective. Note that some have mono- and two-species mixtures, and others could have for instance from 1 to 16 species in their experiments. Considering mono-species litter as a control is therefore the only way to be consistent across all studies. That is, quantification based on a species loss perspective requires us to obtain the effect sizes using mixed-species litter as a control. In this case, different studies have different levels of richness for a control, making it impossible to have a quantification in a standardized way. We first calculated the effect sizes based on mass loss and the decomposition rate constant k (Fig. 1b, c); because of the limited data availability for the constant k, we only calculated the effect sizes for the entire dataset and the subset of data for leaf litter. For mass loss, we further calculated the effect sizes for the different ecosystem types in different climatic regions; subsets of the data that originated from at least three different studies were used. We conducted a multilevel mixed effects meta-regression by considering random effects due to non-independence among some data points, resulting from a nested structure of the dataset that individual data points were nested within a treatment and treatments were also nested within a study. We used the Q statistic for the test of significance.
Note that, because publication bias is a problem in meta-analysis, we visually evaluated the possibility of such a potential bias by plotting the values of the effect sizes and their variances against sample size. If there is no publication bias, studies with small sample sizes should have an increased sampling error relative to those with larger sample sizes, and the variance should decrease with increasing sample size. In addition, the effect size should be independent of the sample size. Also, there should be large variation in effect sizes at the smallest sample sizes. Prior to the meta-analysis, we confirmed that these conditions existed (Supplementary Fig. 1), and found no publication bias large enough to invalidate our analysis.
We also visually confirmed that the datasets were normally distributed based on normal quantile plots. Furthermore, following the method of Gibson et al.70, we randomly selected only one comparison per treatment and then calculated the effect sizes for the decomposition rate (based on mass loss); we repeated this procedure 10,000 times (with replacement) and found that the species-mixing effects on decomposition were significantly positive (0.247 ± 0.045 for the mean ± 95% CI; Supplementary Fig. 2). We thus confirmed that the overall results were not affected by a publication bias or by non-independence of the dataset. We additionally focused on a subset of the entire dataset, which reported mass loss data for each component species from mixtures and thus had comparisons between mass loss of mono-species litter and that of the same species in a mixture (e.g., mass loss of mono-species litter for species A, B, and C was only compared with that for species A, B, and C within a three-species mixture, respectively). We found the results for this subset of the data (Supplementary Fig. 3) almost identical to the main results (Fig. 2a, b). Note that because of the limited data availability, we did not consider the effects of using different equation forms to estimate the rate constant k and instead considered these among-study differences based on a random term in our multilevel meta-analysis. We also calculated the influence of incubation period (days) on the effect size for the decomposition rate (based on mass loss of leaf litter) for the subsets of the dataset using a mixed effects meta-regression30 (Supplementary Fig. 4). We limited this analysis to studies that retrieved litter bags at least two different time points. To account for non-independence of data points, the study identity was included as a random term.
Impacts of global change drivers on litter decomposition
Litter decomposition is primarily controlled by three factors: the environment (e.g., climate), the community of decomposers, and litter traits46,47,48. By relying on the dataset of Makkonen et al.50, who conducted a full reciprocal litter transplant experiment with 16 plant species that varied in their traits and origins (four forest sites from subarctic to tropics), we aimed to obtain a standardized equation to estimate how decomposition rate changed along a climatic gradient (hereafter, the “standardized climate–decomposition relationship”). First, by considering litter decomposition rates at the coldest location (i.e., the subarctic site) as a reference (control), we calculated the effect sizes (i.e., standardized mean difference based on Hedges’ d) for their decomposition rate constant k for all possible comparisons (i.e., between data from the coldest site and the value, one at a time, for the three other warmer sites). As explained above, the effect sizes were calculated only for the comparisons between decomposition rates of litter that were obtained using the same protocol (litter species, origin, and decomposers). The effect sizes therefore represent the climatic effects on litter decomposition after removing the effects of the decomposer community and litter traits.
We obtained the bioclimatic variables for these four study sites from the WorldClim database, and then modelled the standardized climate–decomposition relationships using a multilevel mixed effects meta-regression30,69. Annual mean temperature, the mean temperature of the wettest quarter, or precipitation of the driest quarter were used as an explanatory variable, with the protocol as a random term. This allowed us to evaluate how decomposition rate can be altered by climate along a latitudinal gradient from tropics to subarctic. We then calculated the decomposition rate constant k for our dataset in exactly the same manner as Makkonen et al.50 and calculated the effect sizes; because of limited datasets for most biomes, we performed this analysis only for the subset of our dataset that we obtained for forest biomes, and we used only leaf litter (57 studies) so the results would be comparable with those of Makkonen et al.50. We used these effect sizes and the standardized climate–decomposition relationships to convert the species-mixing effects into climatic effects; that is, based on the temperature or precipitation required to alter the effect size to the same magnitude as the litter diversity effect, we estimated the climate-equivalency of the diversity effects. This means that we relied on the slope of the relationships to quantify the diversity effects, as has been done in previous analyses38,71. We also compared these estimates to the future projections of changes in decomposition rates in response to climate change (with estimates for the 2070s for two of the representative concentration pathway scenarios of CO2 used in the 5th Climate Model Intercomparison Project; CMIP5 RCP 2.6 and 8.5) for the locations of the original studies (the coordinates of individual studies were collected using the Google Map; www.google.com). This comparison was aimed at quantifying the impacts of different global change drivers (biodiversity and climate change) on decomposition processes. Note that we additionally conducted the above analysis after excluding data with no species identity information for mono and/or mixed-species litter bags, those that measured ash-free dry mass loss (this exclusion was to be consistent with the procedure used in Makkonen et al.50), or both of them (Table S1). As an additional confirmation, we also relied on an alternative method used by Hooper et al.45 to convert the species-mixing effects into climatic effects (mean annual temperature), and compared the impacts of different global change drivers on the decomposition rate (Supplementary Fig. 5).
We further analyzed the potential changes in decomposition resulting from litter diversity and climate change. For the dataset of Makkonen et al.50, we first modelled the relationship between climate (annual mean temperature) and the decomposition rate constant k using an LMM with the protocol as a random term. The constant k was log-transformed to improve homoscedasticity. With this equation and the values of annual mean temperature at the 57 study locations in the forest biomes included in our dataset, we then estimated the expected values of k at these study locations (the present k values). Note that because the experiment of Makkonen et al.50 had no mixed-species litter, the expected k using the above LMM was used primarily to estimate the decomposition rate of mono-species litter at a given annual mean temperature. By combining this estimate with the aforementioned estimate of the climate-equivalency of the litter diversity effect (for annual mean temperature; Fig. 3a), it was possible to estimate the potential changes in the k that resulted from increasing litter diversity from mono- to mixed-species litter at the 57 forest study locations. Specifically, the diversity effect converted into the temperature effect (Fig. 3a) was added to present estimates of annual mean temperature at these study locations, and these values of temperature change were used as an explanatory variable to project changes in k based on the above LMM (the projected k values). In other words, we quantified the projected k values under a scenario in which we brought mono-species litter to areas with a warmer annual mean temperature equivalent to the temperature increase in the climate-equivalency effect. We then converted the present and projected values of the decomposition rate constant k to a mass loss per unit time, making it possible to calculate percentage changes in the decomposition rate due to litter diversification. Furthermore, using the above LMM for the temperature–decomposition relationship with the projected changes in mean annual temperature (based on the CMIP5 RCP 2.6 and 8.5 scenarios) as an explanatory variable, we projected future changes in the decomposition rate constant k at these study locations.
Again, we converted these estimates to a litter mass loss to obtain the percentage changes in the decomposition rate that would result from climate change. The rationale for using the temperature–decomposition relationship was as follows: A similar magnitude of increase in k can, in absolute terms, have different influences on decomposition rates at different locations with a different k. For instance, a change in the value of k from 0.1 to 0.2 and a change from 3.0 to 3.1 are substantially different when considered based on litter mass loss or litter amount remaining after decomposition. We thus included annual mean temperature in our model for quantifying the effects of litter diversity on decomposition.