Statistical analysis and structural equation models
All statistical analyses were conducted in R version 4.0.3
(R
Core Team 2020). We calculated the Chao 1 index for the extrapolated
fungal richness (hereafter referred to as fungal richness) of ASVs in
the fungal communities by plot with the estimateR function in thevegan package
(Oksanenet al. 2016). We quantified soil fungal community structure as a
measure of fungal beta diversity by using the three axes of a non-metric
multi-dimensional scaling (NMDS) ordination of Bray–Curtis
dissimilarities of the ASV community matrix using the metaMDSfunction in the vegan package (stress = 0.206; Fig. S2)
(Laliberté
& Legendre 2010). We identified functional guilds of fungi using the
FUNGuild database (Nguyen et
al. 2016). We summed the number of unique ASVs classified as potential
plant pathogenic and ectomycorrhizal (ECM) fungi to calculate richness
within these functional groups.
We used piecewise structural equation models (SEMs) to quantify the
effects of soil fungi in the relationship between the four exogenous
experimental treatments (i.e., tree species richness, functional
diversity, functional identity, and water availability) and annual
above-ground productivity using the piecewiseSEM package
(Lefcheck
2016). Above-ground productivity was square root transformed. The first
set of SEMs were used to quantify the effect of soil fungi on tree
productivity. We developed three models for different fungal components,
i.e., fungal richness, fungal community structure (NMDS axes 1-3), and
functional groups, as well as a null model with no path(s) between
fungal components and productivity. The initial models included the
direct paths of each of the exogenous variables to soil fungal
components and to tree productivity and the direct path(s) from the
fungal components to tree productivity. Indirect effects were assumed to
be significant if their component paths were significant. The SEMs were
constructed with linear mixed effects models with tree species richness,
functional diversity, functional identity, and water treatment as fixed
effects and block as a random effect.
The same approach was used for quantifying effects of the exogenous
treatments on NE and CE using only data for the 128 mixed community
plots. One outlier was removed for the NE model and two from the CE
model based on Grubb’s test (Fig. S3). In all SEMs, we started with the
fully specified model and performed a stepwise elimination of the path
with the highest P-value until the Akaike Information Criterion (AIC)
did not decrease.