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.