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).