Methods

Data collection

Using Google Scholar and Web of Science, we found the peer-reviewed papers through the combination of several search keywords, including (plant diversity OR species diversity OR plant mixture OR species mixture OR mix plant OR polyculture OR intercrop) and (fine root OR root biomass OR root density OR root length density OR root/shoot OR biomass allocation OR specific root length OR SRL OR root diameter OR root nitrogen), up to 1st July 2019. The following criteria were applied for the selection of publications: (1) studies were purposely implemented to isolate the effects of plant species diversity from other factors, such as water treatment and the nutrition addition; (2) values of fine-root traits could be extracted directly from the text, tables, and figures; (3) papers that focused merely on the effects of diversity on root biomass were excluded; (4) genotype mixtures with species were not included; (5) each plant species mixture was compared to corresponding monocultures.
For each study, we extracted the fine-root biomass at different soil depths to calculate the vertical distribution. Fine-root traits were also collected including R/S, RLD, SRL, RD, and RN. For studies that reported root attributes by root order or diameter class, we calculated the community-level means of these values. The plant species richness in mixtures, stand age for forests, or experimental age for grasslands and croplands, and species proportions in mixtures were recorded from the original publications. If different locations, mixture ratios, or abiotic treatments with independent controls were involved in a given publication, we treated them as distinct comparisons (studies) in that publication. In total, 103 studies with 852 paired observations from 64 publications were selected for this meta-analysis.
The proportions of each species in mixtures were based on basal areas or stem densities in forests, the seeds sown in grasslands and croplands, and the number of individuals in containers. Forest stand ages were recorded from the site descriptions in the publications, whereas the experimental ages in containers, grasslands, or croplands were determined by the period between the initiation of the experiment and sampling of the fine roots. Soil sampling depth intervals were converted to the middle values of corresponding depth intervals to facilitate analysis across studies that involved a wide range of depth intervals (Chen & Brassard 2013).
Ecosystem types were categorized as either container, natural forest, planted forest, grassland, or cropland. We obtained geographical locations (altitude, latitude, and longitude) from the original papers that described experiments being conducted in croplands, grasslands, planted forests and natural forests. We recorded the mean annual temperature (MAT) and precipitation (MAP) (when available) conveyed in the original publications or derived based on the geographical location for each site from the WorldClim version 2 Dataset (Fick & Hijmans 2017).

Data analysis

We calculated the community weighted-mean rooting depth (WRD) to compare the fine-root vertical distribution in mixtures with monocultures. WRD was calculated as:
\(\text{WRD\ }\left(\text{cm}\right)=\sum_{i=1}^{n}\left(\frac{B_{i}}{B_{T}}\times D_{i}\right)\)(1)
where Bi is the fine-root biomass in the ith soil layer, BT is the total biomass in all soil layers, and Di is soil sampling depth (as the middle value of each sampling depth interval) of the ith layer.
Natural log-transformed response ratio (lnRR) (Hedges et al. 1999) was employed as the effect size for fine-root biomass and traits (root attributes hereafter). We calculated lnRR as:
\(lnRR=ln\left(\frac{X_{t}}{X_{c}}\right)\) (2)
where Xt is the observed value in the mixture, and Xc is expected value. As Loreau and Hector (2001) recommended, the expected value Xc was calculated as the weighted mean of the corresponding species in monocultures according to the species proportion in mixtures for all root attributes. For root biomass and RLD, Xt is the sum of each constituent species in mixtures. Since R/S, WRD, SRL, RD, and RN are not judged by soil area or volume, Xt was the weighted mean of each constituent species based on the species proportion in mixtures for these traits. For three of the 64 publications in which species proportions were unavailable, we assumed that the species in mixtures were equally distributed. Analysis without the data from these three publications yielded quantitatively similar results. For simplicity and inclusivity, we reported the data from all 64 publications.
As relates to the effect size estimates, we employed the number of replications for weighting (Ma & Chen 2016).
\(W_{r}=\frac{\left(N_{c}\times N_{t}\right)}{\left(N_{c}+N_{t}\right)}\)(3)
where Wr is the weight for each observation, whereas Nt, and Nc are the numbers of replications in mixtures and monocultures, respectively.
To ensure the assumption of linearity between each trait and the species richness in mixtures (R), stand, or experimental age (A), and soil depth (D), we compared the linear, log-linear and quadratic functions for R, A, and D for each root attribute, using equation (4):
\(lnRR=\beta_{0}+\beta_{1}\times X+\pi_{\text{study}}+\varepsilon\)(4)
where β is the estimated coefficient, πstudy is the random effect factor of study, ε is the sampling error, and X is the linear, log-linear, or quadratic form of R, A, and D. We conducted our analysis using the restricted maximum likelihood estimation with the lme4 package with Wr as the weight for each corresponding observation (Bates et al. 2015).
To test the simultaneous effects of R, A, and D on lnRR of each root attribute, we employed the following model:
\(lnRR=\beta_{0}+\beta_{1}\times R+\beta_{2}\times A+\beta_{3}\times D+\beta_{4}\times R\times A+\beta_{5}\times R\times lnD+\beta_{6}\times A\times D+\pi_{\text{study}}+\varepsilon\)(5)
where β is the coefficient to be estimated, πstudy is the random effect factor of study, and ε is the sampling error. The function forms (linear, log-linear, and quadratic) of the three predictors in equation 5 were selected based on the lowest AIC values derived from equation 4 for each root attribute (Table S2). We employed the restricted maximum likelihood estimation with the lme4 package with Wr as the weight for each corresponding observation (Bates et al. 2015). The term D in equation (5) was excluded for R/S and WRD since they are trait variables for the entire ecosystem and entire soil profile, respectively. To prevent overfitting, we derived the most parsimonious model, which was selected using the ‘dredge’ function of the MuMln package (Bartoń 2019).
To elucidate whether the species mixture effects changed with the background environment or ecosystem type, we conducted two types of analysis. First, for those studies conducted in natural habitats, we examined whether lnRR was dependent on MAT and MAP. However, these MAT and MAP effects might be confounded with variations in R, A, and D. Second, we added MAT and MAP and their interactions with the predictors to the most parsimonious models derived from equation 5, from which we then obtained the most parsimonious models to determine whether MAT and MAP accounted for additional variations in lnRR. We also substituted MAT and MAP by ecosystem type and then conducted the same analysis as described above. Moreover, we tested whether the effects of plant mixtures on root attributes associated with MAT and MAP differed between ecosystem types through the linear mixed effect model with ‘study’ as the random effect.
All continuous predictors including R, A, D, MAP, and MAT were scaled to ease the comparisons for fine-root attributes that had variable R, ln(A), ln(D), MAT, and MAP (observed values minus mean and divided by one standard deviation (Cohen et al. 2014). In this way, β0 is the overall mean lnRR at the mean R, mean A, mean D, mean MAP, and mean MAT for each root attribute.
To facilitate interpretation, lnRR and its corresponding 95% confidence interval was transformed to a percentage change using the equation:
\(\left(e^{\text{lnRR}}-1\right)\times 100\%\) (6)
If the CIs did not cover zero, the mixture effect was significant atα  = 0.05. Histograms of model residues and the Shapiro-Wilk test were employed to check the normality of all models, bootstrapped estimates were derived when the normality was violated by using theboot package (Davison & Hinkley 1997; Canty & Ripley 2012). All analyses were performed in R 3.6.1 (R Core Team 2019).