Material and methods

Data collection

We searched for all peer-reviewed publications that investigated the effects of plant diversity on carbon and nitrogen attributes in grasslands using the ISI Web of Science (isiknowledge.com) and China National Knowledge Infrastructure (CNKI, www.cnki.net) databases from 1980-2019. Keyword combinations such as diversity/mixture, grassland/pasture/meadow, carbon/C, nitrogen/N, respiration, and richness were used for the search. The following criteria were applied to select studies: 1) experiments had at least one pair of data comparing monoculture vs. mixture, 2) the species in monoculture were included in the mixtures at the same temporal and spatial scale, 3) response variables (means and measures of variance) were reported for at least one carbon or nitrogen ecosystem function: aboveground biomass (AGB), belowground biomass (BGB), total biomass (TB), soil carbon pool (SCP, %), soil respiration (Rs), heterotrophic respiration (Rh), microbial biomass (MB), bacterial biomass (BB), fungal biomass (FB), aboveground nitrogen pool (ANP), soil nitrogen pool (SNP, %), soil ammonium nitrogen (SAN), soil nitrate nitrogen (SNN), soil nitrogen leaching (SNL), soil nitrogen mineralization (SNM), and 4) the diversity level (species richness), experimental duration, and experimental type (field or greenhouse) were clearly described. Means, standard deviations/errors and samples size were extracted when necessary from digitized graphs using GetData Graph Digitizer (ver. 2.24, < www.getdata-graph-digitizer.com/ >). Sample sizes corresponding to each observation was derived based on the number of independent experimental units. The data were extracted from the Worldclim database at http://worldclim.org/ using the location information (e.g. latitude and longitude). Experimental age was divided into two durations: 1 (1– 4 years) and 4 (>4 years). In total, 73 studies, reporting 15 attributes, and 1385 observations about the effects of plant diversity on carbon and nitrogen processes were selected (Table S1).

Data analysis

We followed the methods of Hedges et al. (1999) and Gurevitch et al. (2018) to evaluate the changes in carbon and nitrogen attributes in plant mixtures. The log response ratio (lnRR , natural log of the ratio of the mean value of monocultures plots to that in mixtures) was calculated as below:
\(lnRR=\ln\left(\frac{{\overset{\overline{}}{x}}_{t}}{{\overset{\overline{}}{x}}_{c}}\right)=\ln\left({\overset{\overline{}}{x}}_{t}\right)-ln({x\overline{}}_{c})\)(1)
where x̅t and x̅c are means of the variable in all mixtures and monocultures plots respectively. Random-effects models were used to estimate weighted average response ratios via the rma function in the metafor package (Cooperet al. 2009) in R 3.5.1 (R Core Team, 2018). Weights for studies were estimates by sampling variance (Hedges et al. 1999) and the between-sampling variability:
\(w_{v}={(\frac{s_{t}^{2}}{n_{t}{\overset{\overline{}}{x}}_{t}^{2}}+\frac{s_{c}^{2}}{n_{c}{\overset{\overline{}}{x}}_{c}^{2}}+\tau^{2})}^{-1}\)(2)
where St, nt, Sc and nc are the standard deviation and sample size for the mixtures and monocultures respectively, and \(\tau^{2}\) the total amount of heterogeneity.
For each attribute, we tested the impacts of species richness (SR), experimental age (EA), and climate using the following model:
lnRR = β0 + β1 × lnSR+ β2 × EA +β3 × climate +β4 × lnSR ×EA+β5 × lnSR ×climate +ɛ (3)
where β is the coefficient to be estimated; ɛ is sampling error. We conducted the analyses using restricted maximum likelihood estimation by the rma function with \(w_{v}\) as the weight for each corresponding observation (Johnson and Omland 2004; Bates et al.2015). For ease of interpretation, lnRR and the corresponding 95% confidence intervals (CIs) were transformed to a percentage change between monocultures and mixtures as (e lnRR-1) ×100%. To illustrate the effects of plant diversity on soil carbon and nitrogen pool with time, we compared the lnRR when the plant diversity in mixtures was R1 (all species present) and Rα (α% lower species richness) using the following equation (Chen et al . 2019):
\(P_{\alpha}={{(R}_{1}/R_{\alpha})}^{\beta_{1}+\beta_{4}T}\) (4)
where Pα is the proportion of remaining soil carbon and nitrogen pool under α% lower species richness in a period ofT. Other model terms were as described for Equation (3). Based on Equation (4), we fitted curves to estimate the potential decrease in soil carbon and nitrogen pool over time under scenarios of 10, 20, 40, and 80% decreases in species richness. We limited the forecasting to 30 years as that is the age of the longest field biodiversity experiment at Cedar Creek Natural History Area, Minnesota, USA (Tilman et al.2006).