Plant mixture balances terrestrial ecosystem C:N:P stoichiometry

Data collection

We systematically searched all peer-reviewed publications that were published prior to May 2021, which investigated the effects of plant diversity on terrestrial C:N:P ratios (i.e., plants, soils, soil microbial biomass, and extracellular enzymes) using the Web of Science (Core Collection; http://www.webofknowledge.com), Google Scholar (http://scholar.google.com), and the China National Knowledge Infrastructure (CNKI; https://www.cnki.net) using the search term: “C:N or C:P or N:P or C:N:P AND plant OR soil OR microbial biomass OR extracellular enzyme OR exoenzyme AND plant diversity OR richness OR mixture OR pure OR polyculture OR monoculture OR overyielding”, and also searched for references within these papers. Our survey also included studies summarized in previously published diversity-ecosystem functioning meta-analyses^15,17,20,33. The literature search was performed following the guidelines of PRISMA (Preferred Reporting Items for Systematic Reviews and Meta-Analyses) (Moher, Liberati⁴⁴; Supplementary Fig. 5).

We employed the following criteria to select the studies: (i) they were purposely designed to test the effects of plant diversity on C:N:P ratios, (ii) they had at least one species mixture treatment and corresponding monocultures, (iii) they had the same initial climatic and soil properties in the monoculture and mixture treatment plots. In thirteen publications, several experiments, each with independent controls, were conducted at different locations and were considered to be distinct studies. In total, 169 studies met these criteria (Supplementary Fig. 5 and Supplementary Table 3). When different publications included the same data, we recorded the data only once. When a study included plant species mixtures of different numbers of species, we considered them as distinct observations.

For each site, we extracted the means, the number of replications, and standard deviations of the C:N, N:P, and C:P ratios of plants (including leaves, shoots, fine roots, total roots), soils, soil enzymes as well as soil microbial biomass C:N ratios, if reported. Similar to Zhou and Staver⁴⁵, we collected nine types of soil enzymes and integrated individual soil enzymes into combined enzymes to represent proxies targeting specific resource acquisitions: C-acquisition (average of Invertase, α-Glucosidase, β-1,4-Glucosidase, Cellobiohydrolase, β-1,4-Xylosidase), N-acquisition (average of β-1,4-N-acetylglucosaminidase, Leucine-aminopeptidase, Urease), and P-acquisition (phosphatase). The ratios of each type of enzyme were subsequently calculated, referred to as soil enzyme C:N, C:P, and N:P. When an original study reported the results graphically, we used Plot Digitizer version 2.0 (Department of Physics at the University of South Alabama, Mobile, AL, USA) to extract data from the figures. This resulted in 52 studies for plant C:N ratios, 35 studies for plant N:P ratios, 17 studies for plant C:P ratios, 83 studies for soil C:N ratios, 42 studies for soil N:P ratios, 19 studies for soil C:P ratios, 33 studies for soil microbial biomass C:N ratios, 41 studies for soil enzyme C:N ratios, 40 studies for soil enzyme N:P ratios and 34 studies for soil enzyme C:P ratios (Supplementary Table 3).

We also extracted species compositions in mixtures, latitude, longitude, stand age, ecosystem type (i.e, forest, grassland, cropland, pot), mean annual temperature (MAT, °C), management practice (fertilization or not), soil type (FAO classification) and sampled soil depth from original or cited papers, or cited data sources. The mean annual aridity index and solar radiation data were retrieved from the CGIAR-CSI Global Aridity Index data set⁴⁶ and WorldClim Version 2⁴⁷ using location information. The annual aridity index was calculated as the ratio of the mean annual precipitation to mean annual potential evapotranspiration⁴⁸. Stand age (SA) was recorded as the number of years since stand establishment following stand-replacing disturbances in forests, and the number of years between the initiation and measurements of the experiments in grasslands, croplands, and pots. Observations were averaged if multiple measurements were conducted during different seasons within a year. The species proportions in plant mixtures were based on the stem density in forests and pots, coverage in croplands, and sown seeds in grasslands. Soil depth was recorded as the midpoint of each soil depth interval⁴⁹. We employed the weighted averages of soil C:N, C:P, and N:P ratios of monocultures in each study as proxies for the status of background nutrients. For studies that did not report soil C:N, C:P, and N:P ratios of monocultures, we used the initial soil C:N, C:P, and N:P ratios (before experiment establishment, if reported) as proxies for the status of background nutrients. When a study reported the soils, soil microbial biomass or soil enzyme C:N:P data from multiple soil depths, we used the soil C:N, C:P, and N:P ratios of the corresponding depths as background nutrient proxies. For plant C:N:P data, we used the uppermost soil layer C:N, C:P, and N:P ratios as background nutrient proxies, since it contains the majority of the available nutrients essential for plant growth⁵⁰. We compared the estimates for the data sets with and without pot studies and found that both data sets yielded qualitatively similar results (Supplementary Tables 2 and 4). Thus, we reported results based on the whole data set.

We employed two key functional traits to describe the functional composition: ‘leaf nitrogen content per leaf dry mass’ (N_mass, mg g⁻¹), and “specific leaf area” (SLA, mm² mg⁻¹; i.e., leaf area per leaf dry mass), as they are expected to be related to plant growth rate, resource uptake and use efficiency²⁷, and are available for large numbers of species. We obtained the mean trait values of N_mass and SLA data by using all available measurements for each plant species from the TRY Plant Trait Database⁵¹ except for two studies that included the data in their original publication⁵², or related publications in the same sites⁵³. Functional diversity (FD_is) was calculated as functional dispersion, which is the mean distance of each species to the centroid of all species in the functional trait space, based on the two traits together⁵⁴. The calculation of FD_is was conducted using the FD package⁵⁴.

Data analysis

The natural log-transformed response ratio (lnRR) was employed to quantify the effects of plant mixture following Hedges, Gurevitch⁵⁵:

where ({overline{X}}_{{{{{{rm{t}}}}}}}) and ({overline{X}}_{{{{{{rm{c}}}}}}}) are the observed values of a selected variable in the mixture and the expected value of the mixture in each study, respectively. If a study has multiple richness levels in mixtures (for example, 1, 4, 8, and 16), lnRR was calculated for the species richness levels 4, 8, and 16, respectively. We calculated ({overline{X}}_{{{{{{rm{c}}}}}}}) based on weighted values of the component species in monocultures following Loreau and Hector³⁹:

where m_i is the observed value of the selected variable of the monoculture of species i and p_i is the proportion of species i density in the corresponding mixture. When a study reported multiple types of mixtures (species richness levels) and experimental years, ({overline{X}}_{{{{{{rm{t}}}}}}}) and ({overline{X}}_{{{{{{rm{c}}}}}}}) were calculated separately for each mixture type and experimental year.

In our data set, sampling variances were not reported in 37 of the 169 studies, and no single control group mean estimate is present with standard deviation or the standard error reported. Like the previous studies^6,56, we employed the number of replications for weighting:

where W_r is the weight associated with each lnRR observation, and N_c and N_t are the number of replications in monocultures and corresponding mixtures, respectively.

The C:N, N:P, and C:P ratios of plants, soils, and soil enzymes, as well as soil microbial biomass C:N ratios were considered as response variables and analyzed separately. To validate the linearity assumption for the continuous predictors, we initially graphically plotted the lnRR vs. individual predictors, including FD_is, SA, and background nutrient status (N, i.e., C:N, C:P, and N:P ratios of soil) and identified logarithmic functions as an alternative to linear functions. We also statistically compared the linear and logarithmic functions with the predictor of interest as the fixed effect, and “study” and measured plant parts (i.e., leaves, shoots, fine roots, total roots) or soil depth as the random effects, using Akaike information criterion (AIC). The random factors were used to account for the autocorrelation among observations within each “Study”, and potential influences of variation in measured plant parts and soil depth. We found that the linear FD_is, SA, and N resulted in lower, or similar AIC values (∆AIC < 2) except for plant C:N and enzyme C:P (Supplementary Table 5). The logarithmic N and FD_is resulted in lower AIC values for plant C:N and soil enzyme C:P, respectively (Supplementary Table 5). We used the following Eq. 4 to determine the effects of the FD_is, SA, and environmental variables (E, i.e., background nutrient status, ecosystem type, the proportion of N-fixing plants, mean annual solar radiation (S), and aridity index (AI)) and their interactions on the C:N, N:P, C:P ratios of plants, soils, soil microbial biomass, and enzymes:

where β_i and ɛ are coefficients and sampling error, respectively; π_study was random effects accounting for the autocorrelation among observations within each “Study” while π_p and π_depth were random effects accounting for the potential influences of variation in measured plant parts (i.e., leaves, shoots, fine roots, total roots, only used for plant C:N:P) and sampled soil depth (used for soil, soil microbial biomass, and enzyme C:N:P as well as plant C:N:P when N was included as a predictor), respectively. We conducted the analysis using the restricted maximum likelihood estimation with the lme4⁵⁷ and lmerTest⁵⁸ package. We centered all continuous predictors (observed values minus mean). When continuous predictors were scaled, β₀ is the overall mean lnRR at the mean FD_is and SA and N⁵⁹.

To prevent overfitting⁶⁰, we selected the most parsimonious model among all alternatives instead of using stepwise multiple regression, which can be biased and has multiple shortcomings⁶¹. We employed the condition of retaining the FD_is as it was intrinsic to the purpose of the study for assessing the effects of plant diversity in mixtures. Model selection among alternatives was accomplished using the “dredge” function of the muMIn package⁶². All terms associated with stand age and ecosystem type, the proportion of N-fixing plants or S were excluded in the most parsimonious models. The model selection led to Eq. 5 for the plant C:N, C:P and soil C:N, N:P, C:P ratios, Eq. 6 for the plant N:P, and enzyme C:N, N:P, C:P ratios, and Eq. 7 for the microbial biomass C:N, respectively:

As an alternative to functional diversity, we also investigated how plant and soil C:N:P changed with species richness, by replacing FD_is in Eqs. 5, 6, and 7 with species richness. We found that both diversity metrics yielded qualitatively similar trends for all C:N:P variables, except soil C:N and C:P, which significantly decreased with species richness (Supplementary Fig. 2). However, these results are driven by a single study (Jena) with high species richness (4, 8, 16) or limited species richness levels for some selected variables (2 levels soil C:P), respectively (Supplementary Fig. 2).

To further examine the effects of environmental variables, we conducted an analysis with the environment variable (ecosystem type, the proportion of N-fixing plants, mean annual solar radiation, aridity index, soil types, background soil nutrient status, and management practice) as the only fixed factor, and ‘study’ and measured plant parts as the random factors. The analysis confirmed that there were no differences in the responses of the C:N, N:P, and C:P ratios of plants and soils, as well as soil microbial biomass C:N ratios to mixtures among experimental environment variables (Supplementary Table 1).

We employed partial regressions (or predicted effects) to graphically demonstrate the effects of predictors. Our analysis indicated that many of our models violated the assumption of normality, based on Shapiro–Wilk’s test on model residuals. Consequently, we bootstrapped the fitted coefficients by 1000 iterations⁶³ with the boot package⁶⁴. The collinearity among explanatory variables was examined by evaluating the variance inflation factor (VIF, only models with all predictors having VIFs < 3 were accepted⁶⁵), and no multicollinearity problem was found in the most parsimonious models (Supplementary Table 2). We analyzed the potential for publication bias to influence our results using Egger’s regressions test for funnel plot asymmetry on mixed-effects models⁶⁶, with sample size as the predictor. Egger’s test was run on the main statistical tests we performed (response ratio across the entire data set, and then the response ratio including associated predictors as covariates). We did not find any significant publication bias that could bias our results toward significant effects according to Egger’s regression (Supplementary Table 6).

The coefficients or treatment effects were significant from zero at α = 0.05 if the bootstrapped 95% confidence intervals (CIs) did not cover zero. The mean effect sizes between groups were significantly different if their 95% CIs did not overlap the other’s mean. To facilitate interpretation, we transformed the lnRR and its corresponding CIs back to a percentage using (({e}^{{{{{{{mathrm{ln}}}}}}}{RR}}-1)) × 100%. Furthermore, we handled data using data.table⁶⁷, and visualized data using maps⁶⁸ and ggplot2⁶⁹. All package versions are provided in the Reporting summary. All statistical analyses were performed in R 4.0.2⁷⁰.

Reporting summary

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

Source: Ecology - nature.com