Study sites and sampling design<\/p>\n
We used a pan-European network of 64 mature, uneven-aged forest plots (30\u2009\u00d7\u200930\u2009m2) consisting of three-species mixture stands (34 plots) and corresponding monospecific stands (30 plots; Extended Data Table 1<\/a>). These plots are part of the FunDivEUROPE\u00a0(Functional Significance of Forest Biodiversity in Europe) exploratory platform59<\/a>, and were established across European forests over 2011\u20132012 to investigate the role of the diversity and composition of regionally common and economically important tree species on ecosystem functioning. The studied plots were distributed across four locations featuring different European forest types and spanning a large biogeographic gradient: North Karelia (Finland), Bia\u0142owie\u017ca (Poland), R\u00e2\u015fca (Romania) and Colline Metallifere (Italy), corresponding to typical boreal, hemi-boreal, mountainous beech and thermophilous deciduous (Mediterranean) forests, respectively (Supplementary Table 4<\/a> and Supplementary Fig. 5a<\/a>). In each location, plots were carefully selected based on tree species richness and composition while minimizing as much as possible covariation with potentially confounding environmental factors such as topography and soil conditions59<\/a> (Supplementary Fig. 2a<\/a>). Plot selection was performed so as to include monospecific stands of all tree species from the local species pool and replicate the three-species mixture treatment with different tree species combinations while maximizing community evenness (Extended Data Table 1<\/a>). This allowed strict avoidance of a dilution gradient, such as would occur in a design with monospecific stands of only one species combined with mixture stands including this species, along with a clear distinction between the effects of species mixing and composition. Our stratified plot selection procedure enabled us to mimic formal biodiversity experiments, given that such manipulative approaches are virtually impossible to undertake in mature forests owing to the high longevity of tree species. Tree species diversity and composition in the studied plots were predominantly the result of natural community assembly from the regional species pool, combined with local forest management practices. The investigated levels of species richness, that is, one versus three tree species, are typical for European forest ecosystems (https:\/\/forest.eea.europa.eu\/topics\/forest-biodiversity-and-ecosystems\/forest-ecosystems<\/a>), boreal forests in Asia and North-America, as well as managed forests and plantations worldwide, although clearly much less diversified than many (sub)tropical forests and some mid-latitude temperate forests47<\/a>. Overall, our sampling design encompassed a total pool of 13 tree species, including 12 ectomycorrhizal and one arbuscular mycorrhizal tree species, and the local species pool ranged from three to five tree species per location (Extended Data Table 1<\/a> and Supplementary Table 4<\/a>).<\/p>\n Soil organism sampling and analysis<\/p>\n In each plot, we assessed energy fluxes through the soil food web by measuring the biomass of major groups of organisms within this food web22<\/a>,23<\/a>,43<\/a>,60<\/a>, including both microbial (bacteria and fungi) and faunal (nematodes, microarthropods and macroinvertebrates) groups. Biomass data for all soil organism groups were expressed per unit surface area (g dry weight\u2009m\u22122) at the plot level. Further details on biomass calculation and the trophic classification of soil organisms are provided in the Supplementary Methods<\/a>.<\/p>\n Sampling<\/p>\n Soil organisms were sampled in all plots during the phenological spring of 2017 (Supplementary Table 4<\/a>), a period of high soil biological activity. We sampled both the litter layer (unfragmented aboveground litter, OL horizon) and the soil layer (including both fragmented\/humified organic matter and mineral soil, OF\/OH\/A horizons). In each plot, we selected five 10\u2009\u00d7\u200910-m2 subplots, with samples taken equidistantly from three trees of either the same species in monospecific stands or of different species for three-species mixture stands (Supplementary Fig. 5b<\/a>). For each subplot, soil samples for microbial analyses and nematode extraction were collected by taking five soil cores (10-cm depth, 5.3-cm diameter), spaced approximately 35\u2009cm apart around the equidistant point between the three trees, weighted by tree individual size, that is, individual diameter at breast height. The five cores were gently sieved through 6-mm mesh (to avoid damaging nematodes), homogenized and pooled at the subplot level for nematode extraction. Pooled soil was then sieved through 2-mm mesh for microbial analyses. All subsamples were stored at 4\u2009\u00b0C until further processing. For microarthropod extraction, an intact core (10-cm depth, 10-cm diameter), including both the litter and soil layers, was collected from each of three subplots along a southwest\u2013northeast transect and stored at 4\u2009\u00b0C until further processing. For the hand-sorting of soil macroinvertebrates, an intact monolith (25-cm depth, 25\u2009\u00d7\u200925-cm2 surface), including both the litter and soil layers, was collected from each of the same three subplots. To express all data per unit surface area, an extra core was sampled, sieved through 2-mm mesh, dried at 105\u2009\u00b0C for 48\u2009h and weighted to measure soil bulk density.<\/p>\n Microorganisms<\/p>\n The biomass of bacteria (gram-positive and gram-negative), arbuscular mycorrhizal fungi and non-arbuscular mycorrhizal fungi was quantified at the plot level using phospholipid fatty acid data61<\/a>. Fungal community data based on metagenomic amplicon sequencing and bioinformatics62<\/a> were used to partition non-arbuscular mycorrhizal fungal biomass into five trophic guilds: ericoid mycorrhizal fungi, ectomycorrhizal fungi, general saprotrophic fungi, wood saprotrophic fungi and plant pathogenic fungi (Supplementary Methods<\/a>). The biomass of each fungal trophic guild was calculated by multiplying its relative abundance, that is, number of reads divided by the total number of reads for the five trophic guilds, by the total non-arbuscular mycorrhizal fungal biomass.<\/p>\n Nematodes<\/p>\n Nematodes were extracted for each subplot within 72\u2009h after sampling from approximately 100\u2009g of fresh soil using a modified sugar flotation method63<\/a>, before being heat-killed and fixed in 4% formaldehyde. Nematodes were then pooled and counted at the plot level, and a subsample of approximately 160 randomly selected individuals were identified to family level. The biomass of nematode families was calculated based on body mass data retrieved from the Nemaplex database (http:\/\/nemaplex.ucdavis.edu<\/a>). Nematode families were assigned to five trophic guilds64<\/a>: herbivores, bacterivores, fungivores, omnivores and carnivores.<\/p>\n Microarthropods<\/p>\n Microarthropods were extracted from the intact core (including both the litter and soil layers) within 72\u2009h after sampling for each subplot using the Berlese\u2013Tullgren funnel method65<\/a>, and were fixed in 70% ethanol. Microarthropods were then counted and identified to species level for collembola and to order level for mites. Microarthropod biomass was estimated based on an allometric model using body length data retrieved from the BETSI database (https:\/\/portail.betsi.cnrs.fr<\/a>) for collembola, and data from the literature for mites60<\/a>. Microarthropod taxa were assigned to seven trophic guilds64<\/a> belonging to the following broad trophic groups: microbi-detritivores, fungivores, omnivores and carnivores.<\/p>\n Macroinvertebrates<\/p>\n Macroinvertebrates were hand sorted in the field for each subplot and fixed in 70% ethanol. Macroinvertebrates were then counted and identified to species level for Lumbricidae, Isopoda, Diplopoda, Chilopoda and Araneae; and order to family level for other taxa66<\/a>. All macroinvertebrate individuals were weighed for body mass. Macroinvertebrates were assigned to 23 trophic guilds64<\/a> belonging to the following broad trophic groups: herbivores, detritivores, humi-detritivores, omnivores and carnivores.<\/p>\n Metabolic rates<\/p>\n We used soil microbial respiration and biomass data67<\/a> to calculate the metabolic rates of microbial groups, whereas we used a model based on individual body mass, environmental temperature (mean soil temperature during the growing season; see \u2018Microclimate\u2019 section below) and phylogenetic grouping68<\/a> to calculate the metabolic rates of faunal groups (Supplementary Methods<\/a>).<\/p>\n Assimilation efficiencies<\/p>\n Assimilation efficiencies (that is, the proportion of consumed food assimilated by digestion) specific to each food type (plant-derived resource or prey) for faunal consumers were calculated using a model based on food N content69<\/a>, and we also applied a temperature correction of assimilation efficiency70<\/a> (Supplementary Methods<\/a>). Assimilation efficiencies of all food resources were set to 1 for microbial consumers, given their external digestion system.<\/p>\n Plant-derived resources<\/p>\n In each plot, we quantified the biomass of the following plant-derived (basal) resources of the soil food web: living fine roots38<\/a> (absorptive roots belonging the first three root orders) and associated photosynthates\/rhizodeposits34<\/a>, litter (including dead leaves67<\/a>, dead roots71<\/a> and dead wood) and SOM72<\/a> (Supplementary Methods<\/a>). Photosynthates\/rhizodeposits refer to organic compounds provided by roots to mycorrhizal fungi, or released directly by roots into the soil. To estimate the biomass of rhizodeposits, we used a mass ratio of 0.4 between net rhizodeposition (the portion of rhizodeposited C remaining in the soil after microbial utilization and respiration) and root biomass, based on the results of a meta-analysis of fixed C partitioning in plant\u2013soil systems34<\/a>. The biomass of rhizodeposits was calculated by multiplying the biomass of living fine roots by this factor.<\/p>\n Food web reconstruction<\/p>\n We constructed our soil food web from 47 network nodes, including six plant-derived (basal) resources and 41 trophic guilds of soil organisms (consumers) that were differentiated based on multiple traits64<\/a> (Supplementary Table 5<\/a>). To establish the topology of the soil food web and quantify trophic interaction strengths (Supplementary Fig. 1<\/a>), a food web interaction matrix was constructed based on basic food web principles, and a priori knowledge of soil organism biology and key traits of consumers following the approach in ref.\u20097<\/a>, except that microbes were considered here as consumers rather than basal resources (Supplementary Methods<\/a>). Briefly, the food web interaction matrix was calculated by multiplying five matrices representing different trait dimensions: phylogenetically defined feeding preferences, density dependence, predator\u2013prey interactions related to body mass ratio and hunting strategy, prey defence mechanisms and spatial niche overlap related to vertical stratification. The five matrices relied on the following assumptions, respectively7<\/a>: (1) there are phylogenetically conserved differences in feeding preferences of consumers; (2) food consumption is density (biomass) dependent, that is, consumers will preferentially feed on food resources that are locally abundant owing to a higher encounter rate; (3) the strength of predator\u2013prey interactions is primarily defined by the optimum predator\u2013prey mass ratio, that is, a predator is typically larger than its prey, but certain predator hunting traits can modify the optimum predator\u2013prey mass ratio; (4) the strength of predator\u2013prey interactions can be weakened by prey defence traits, that is, prey with efficient physical, chemical or behavioural protection are consumed less; (5) the strength of trophic interactions between a consumer and its food resource is modulated by the overlap in their spatial niches related to vertical stratification, with greater overlap leading to a stronger interaction. Food web reconstruction was carried out separately for each plot to account for plot-specific density dependence.<\/p>\n Food web energy fluxes and functions<\/p>\n For the calculation of energy fluxes, we assumed a steady state of the soil food web6<\/a>. This means that the energy flowing into a given feeding guild of the food web through food consumption balances the energy lost by excretion, metabolism and predation of that feeding guild. Energy fluxes to each feeding guild within the food web (kJ\u2009m\u22122\u2009d\u22121) were then calculated based on the trophic interaction matrix using the following equation6<\/a>: \\(F=\\frac{1}{{e}_{{\\rm{a}}}}\\times (X+L)\\), where F is the total flux of energy into the feeding guild, ea is the diet-specific assimilation efficiency, X is the community metabolism of the feeding guild and L is the energy loss to predation (see the Supplementary Methods<\/a> for further details). To simplify the representation of the food web, we aggregated biomass and energy flux matrices at broader trophic group levels by summing the rows and columns of trophic nodes belonging to the same trophic group (Fig. 1<\/a> and Supplementary Table 5<\/a>).<\/p>\n We then calculated five broad trophic functions of the soil food web (Fig. 1<\/a>): carnivory, the sum of energy fluxes outgoing from fauna to their faunal consumers; microbivory, the sum of energy fluxes outgoing from microbes to their faunal consumers; herbivory, the sum of energy fluxes outgoing from living fine roots to their consumers; plant C allocation to soil by living roots, the sum of energy fluxes outgoing from photosynthates and rhizodeposits (via living fine roots) to their microbial consumers; and detritivory, the sum of energy fluxes outgoing from detritus (dead organic matter, including plant litter and SOM) to their consumers. Additionally, we calculated eight more specific trophic functions (Fig. 1<\/a>): bacterivory and fungivory, the sum of energy fluxes outgoing from bacteria and fungi to their faunal consumers, respectively; root pathogenicity and rhizophagy, the sum of energy fluxes outgoing from living fine roots to their microbial and faunal consumers, respectively; litter decomposition and litter engineering, the sum of energy fluxes outgoing from plant litter to their microbial and faunal consumers, respectively; and SOM decomposition and soil engineering, the sum of energy fluxes outgoing from SOM to their microbial and faunal consumers, respectively. Decomposition refers to the assimilation and mineralization of dead organic matter through respiration, a process that is mediated mainly by microbes in soil60<\/a>. Engineering refers to the physical modification, maintenance or creation of habitats1<\/a>, which is a major way through which soil faunal detritivory affects the decomposition, transformation and translocation of dead organic matter7<\/a>,73<\/a>. Such inferences from energy fluxes about effects that are not purely trophic are especially justifiable in soil where habitat and food resources are tightly interlinked7<\/a>. As such, the faunal consumption of litter results in the conversion of litter into faeces which in turn accelerates its decomposition through chemical and physical changes73<\/a>, as a form of litter engineering. Similarly, the faunal consumption of SOM is linked to biopedturbation and soil structure maintenance73<\/a> in a manner that represents soil engineering.<\/p>\n We quantified soil food web multifunctionality based on ten trophic functions: plant C allocation to soil by living roots, root pathogenicity, rhizophagy, litter decomposition, litter engineering, SOM decomposition, soil engineering, bacterivory, fungivory and carnivory. Defining high multifunctionality by fast rates of multiple functions simultaneously, we calculated soil food web multifunctionality using the \u2018averaging\u2019 approach74<\/a>, that is, the main range-standardized values for the ten trophic functions (Supplementary Methods<\/a>). Consistent results were found when using the alternative \u2018threshold\u2019 approach74<\/a> (see \u2018Sensitivity analyses\u2019 in the Supplementary Methods<\/a>). Additionally, we adopted the \u2018single functions\u2019 approach to help illuminate which individual functions drive trends in the effects of tree communities on soil food web multifunctionality74<\/a>.<\/p>\n Tree functional (trait-based) properties<\/p>\n Functional diversity and composition of tree communities were characterized using a set of nine plant functional traits known to be directly involved in resource economics9<\/a>,32<\/a>. These traits included three leaf traits: leaf nitrogen content, specific leaf area and leaf dry matter content; and six fine (absorptive) root traits71<\/a>: root nitrogen content, root tissue density, specific root length, mean root diameter, root length density and ectomycorrhizal colonisation intensity. Trait data of each tree species were mostly derived from plot-specific measurement in the field, but specific leaf area and leaf dry matter content data for Poland and Italy were retrieved from the TRY plant trait database75<\/a> (https:\/\/www.try-db.org<\/a>). Using trait values from databases has limitations because of intraspecific trait variability. However, the sensitivity analysis that we performed revealed that the associated uncertainty in trait values had little bearing on our overall inferences (see \u2018Sensitivity analyses\u2019 in the Supplementary Methods<\/a>). See the Supplementary Methods<\/a> for further details and Supplementary Table 6<\/a> for trait values of each tree species and geographic location.<\/p>\n To quantify functional diversity, we calculated the functional dispersion (FDis) index corresponding to the mean distance of each species to the centroid of all species within the multidimensional trait space76<\/a>. For functional composition, we calculated CWMs of each trait. Both FDis and CWM values were computed based on the relative basal area of tree species. We performed a principal component analysis (PCA) on all CWM traits, which simplified functional composition into two dimensions71<\/a> (Extended Data Fig. 4a<\/a> and Supplementary Table 7<\/a>): (1) an LES (45.2% of variation), which ranged from slow\/conservative leaf attributes (N-poor and dry-matter-rich leaves with low leaf area per unit mass) to fast\/acquisitive leaf attributes (N-rich and dry-matter-poor leaves with high leaf area per unit mass), and also aligned here with a fine-root gradient of belowground resource foraging strategies ranging from low foraging efficiency (thick fine roots with low length per unit mass) to high foraging efficiency (thin fine roots with high length per unit mass); (2) an RES (31.5% of variation), which ranged from slow\/conservative fine-root attributes (N-poor fine roots with low tissue density) to fast\/acquisitive fine-root attributes (N-rich fine roots with high tissue density), and also aligned here with a fine-root gradient of soil exploration strategies ranging from \u2018do-it-yourself\u2019 attributes (high root length density and low ectomycorrhizal colonisation intensity) to \u2018outsourcing\u2019 attributes (low root length density and high ectomycorrhizal colonisation intensity). The mean values of FDis and the scores of the two first PCA axes of CWM trait ordination for each location are shown in Extended Data Table 1<\/a>.<\/p>\n Tree and understorey vegetation<\/p>\n We characterized forest vegetation using a set of eight variables measured in each plot: tree aboveground biomass, tree aboveground litterfall, total root biomass, aboveground biomass of both woody and herbaceous understorey plants, and aboveground C:N ratio and species richness of understorey plant communities (Supplementary Methods<\/a>).<\/p>\n Environmental driversAbiotic conditions<\/p>\n We characterized abiotic conditions in each plot using a set of six variables related to soil texture, macroclimate and topography: the sand, silt and clay content of mineral soil (A horizon, 10-cm depth), mean annual temperature and precipitation, and altitude (Supplementary Methods<\/a>). We performed a PCA to reduce the dimensionality of abiotic properties into two dimensions (Supplementary Fig. 2a<\/a>): (1) a soil texture and macroclimate gradient ranging from coarse soil texture, and dry and cold macroclimate to fine soil texture, and wet and hot macroclimate (70.8% of variation); (2) a temperature gradient ranging from cold macroclimate with high elevation to warm macroclimate with low elevation (19.1% of variation).<\/p>\n Microclimate<\/p>\n We characterized microclimate using a set of six variables: soil temperature and moisture (measured at 8-cm soil depth) and air temperature (measured 50\u2009cm above the ground) for both annual and growing season (daily mean temperatures\u2009>\u20095\u2009\u00b0C) periods (Supplementary Methods<\/a>). We performed a PCA to simplify microclimatic properties into a single dimension ranging from cold and wet to hot and dry (Extended Data Fig. 4b<\/a>; 63.1% of variation).<\/p>\n Leaf litter quality<\/p>\n We characterized freshly fallen tree leaf litter quality using a set of 16 chemical and physical variables mainly related to elemental composition and C quality77<\/a>: the concentrations of N, P, Ca, Mg and K, and the C:N and C:P ratios, the proportion of C fractions (lignin, cellulose, hemicellulose and water-soluble compounds), the concentrations of secondary metabolites (condensed tannins, total phenolics and soluble phenolics), as well as the litter pH and water holding capacity. Leaf litter quality data of each tree species were derived from location-specific measurement (Supplementary Methods<\/a>).<\/p>\n We then quantified the functional diversity and composition of tree leaf litter by calculating the FDis index and CWMs of each litter property based on the relative leaf litter mass of the component tree species. We performed a PCA to simplify tree leaf litter properties into one dimension ranging from low to high nutritional quality (Extended Data Fig. 4c<\/a>; 42.5% of variation).<\/p>\n Soil fertility<\/p>\n We characterized soil fertility using a set of six variables measured in each plot72<\/a>: the organic C content of mineral soil (A horizon, 10-cm depth), the mass of the forest floor (including unfragmented aboveground litter and fragmented\/humified organic matter, OL\/OF\/OH horizons), as well as the pH and C:N ratios of both the forest floor and mineral soil (Supplementary Methods<\/a>). We performed a PCA to simplify soil properties into a single dimension ranging from low to high fertility (Supplementary Fig. 2b<\/a>; 40.1% of variation).<\/p>\n Measured ecosystem processes<\/p>\n To characterize in situ patterns of litter decomposition, we used field-based data of naturally occurring tree leaf litter decomposition measured in each plot78<\/a> (Supplementary Methods<\/a>). To characterize SOM decomposition, we further used soil microbial respiration and biomass data67<\/a> (Supplementary Methods<\/a>). To characterize soil and litter engineering, we used humus type and forest floor mass data72<\/a> as these two variables reflect the transformation and translocation of dead organic matter by faunal detritivory73<\/a>.<\/p>\n Statistical analyses<\/p>\n All analyses were performed using R v.4.3.0 (ref.\u200979<\/a>). To test how tree community properties affect soil food web functioning, we used Bayesian multi-level models that accounted for the hierarchical design of the study, with plots nested within four geographic locations across Europe. We used random intercept models (with varying intercepts but a common slope across locations) to assess general patterns across European forests. Similar results were found when using random slope models (with varying intercepts but a common slope across locations), indicating that our findings are robust to model specification (see \u2018Sensitivity analyses\u2019 in the Supplementary Methods<\/a>). To ensure that the effects of tree communities were not biased by confounding factors, we also explicitly included abiotic covariates into the models, that is, the first two PCA axes representing abiotic conditions (Supplementary Fig. 2a<\/a>). This statistical adjustment allowed us to control for variation in abiotic conditions both among and within locations. Following the principle of the \u2018structural causal model\u2019 framework, controlling for confounding factors allows to satisfy the \u2018backdoor criterion\u2019 required for quantifying unbiased causal relationships from observational data80<\/a>. We also checked whether omitted variable bias generated by potential unmeasured confounders at the location level could affect our inference, and found that it was not an issue in our study (see \u2018Sensitivity analyses\u2019 in the Supplementary Methods<\/a>). We adopted both taxonomic and functional approaches. For the taxonomic approach, we used linear mixed-effects models\u00a0(LMMs) including tree species richness (comparing three-species mixture stands to corresponding monospecific stands) and abiotic covariates (see above) as fixed factors, and tree species composition (a factor listing all tree species present; Extended Data Table 1<\/a>) and location as random factors, as in the following equation:<\/p>\n $$\\begin{array}{c}{y}_{ij}=\\alpha +{\\beta }_{1}{{\\rm{r}}{\\rm{i}}{\\rm{c}}{\\rm{h}}{\\rm{n}}{\\rm{e}}{\\rm{s}}{\\rm{s}}}_{ij}+{\\beta }_{2}{{\\rm{a}}{\\rm{b}}{\\rm{i}}{\\rm{o}}{\\rm{t}}{\\rm{i}}{\\rm{c\\; PC}}1}_{ij}+{\\beta }_{3}{{\\rm{a}}{\\rm{b}}{\\rm{i}}{\\rm{o}}{\\rm{t}}{\\rm{i}}{\\rm{c\\; PC}}2}_{ij}+{\\delta }_{ij}^{{\\rm{c}}{\\rm{o}}{\\rm{m}}{\\rm{p}}{\\rm{o}}{\\rm{s}}{\\rm{i}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}}+{\\delta }_{i}^{{\\rm{l}}{\\rm{o}}{\\rm{c}}{\\rm{a}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}}+{{\\epsilon }}_{ij}\\\\ {\\rm{w}}{\\rm{i}}{\\rm{t}}{\\rm{h}}\\,{\\delta }_{ij}^{{\\rm{c}}{\\rm{o}}{\\rm{m}}{\\rm{p}}{\\rm{o}}{\\rm{s}}{\\rm{i}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}} \\sim N(0,{{\\sigma }}_{{\\rm{c}}{\\rm{o}}{\\rm{m}}{\\rm{p}}{\\rm{o}}{\\rm{s}}{\\rm{i}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}}^{2}),{\\delta }_{i}^{{\\rm{l}}{\\rm{o}}{\\rm{c}}{\\rm{a}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}}\\, \\sim N(0,{{\\sigma }}_{{\\rm{l}}{\\rm{o}}{\\rm{c}}{\\rm{a}}{\\rm{t}}{\\rm{i}}{\\rm{o}}{\\rm{n}}}^{2})\\\\ \\,{\\rm{a}}{\\rm{n}}{\\rm{d}}\\,{{\\epsilon }}_{ij} \\sim N(0,{{\\sigma }}^{2})\\end{array}$$<\/p>\n \n (1)\n <\/p>\n where i refers to location and j refers to plot, y is the dependent variable, \u03b1 is the general intercept, \u03b2x terms are partial slopes, \u03b4 terms are random effects, that is, factor-specific deviation from the common intercept, \u03f5 is the model error, that is, residuals, and \u03c32 is the variance. For the functional approach, we used LMMs including tree functional diversity (FDis), tree functional composition (first two PCA axes of tree CWM traits corresponding, respectively, to the LES and RES; Extended Data Fig. 4a<\/a>) and abiotic covariates (see above) as fixed factors, and location as a random factor, as in the following equation:<\/p>\n $$\\begin{array}{c}{y}_{{ij}}=\\alpha +{\\beta }_{1}{\\mathrm{FDis}}_{{ij}}+{\\beta }_{2}{\\mathrm{CWM\\; trait\\; PC}1}_{{ij}}+{\\beta }_{3}{\\mathrm{CWM\\; trait\\; PC}2}_{{ij}}\\\\ \\,+{\\beta }_{4}{\\mathrm{abiotic\\; PC}1}_{{ij}}+{\\beta }_{5}{\\mathrm{abiotic\\; PC}2}_{{ij}}+{\\delta }_{i}^{\\mathrm{location}}+{\\varepsilon }_{{ij}}\\\\ \\,\\mathrm{with}\\,{\\delta }_{i}^{\\mathrm{location}} \\sim N(0,{\\sigma }_{\\mathrm{location}}^{2})\\,\\mathrm{and}\\,{{\\epsilon }}_{{ij}} \\sim N(0,{\\sigma }^{2})\\end{array}$$<\/p>\n \n (2)\n <\/p>\n To further investigate the effects of functional composition, we also used LMMs including individual tree CWM trait and abiotic covariates as fixed factors, and location as a random factor, as in the following equation:<\/p>\n $$\\begin{array}{c}{y}_{{ij}}={\\beta }_{0}+{\\beta }_{1}{\\mathrm{CWM\\; trait}}_{{ij}}+{\\beta }_{2}{\\mathrm{abiotic\\; PC}1}_{{ij}}\\\\ \\,+{\\beta }_{3}{\\mathrm{abiotic\\; PC}2}_{{ij}}+{\\delta }_{i}^{\\mathrm{location}}+{\\varepsilon }_{{ij}}\\\\ \\mathrm{with}\\,{\\delta }_{i}^{\\mathrm{location}} \\sim N(0,{\\sigma }_{\\mathrm{location}}^{2})\\,\\mathrm{and}\\,{{\\epsilon }}_{{ij}} \\sim N(0,{\\sigma }^{2})\\,\\end{array}$$<\/p>\n \n (3)\n <\/p>\n Bayesian LMMs were fitted with Markov chain Monte Carlo methods with non-informative priors for all parameters. Each model was fitted based on 10,000 iterations of both warm-up and sampling phases and was checked for their convergence and compliance with statistical assumptions, including normality and independence of residuals, normality of random effects, homogeneity of variance, linearity of relationships and absence of multicollinearity (Supplementary Methods<\/a>).<\/p>\n For each model, we performed variance partitioning following a model-based approach based on variance components81<\/a> to quantify the proportion of variance in the dependent variable explained by three groups of predictors: diversity (tree species richness or FDis), composition (tree species composition or CWM traits PC1 and PC2) and biogeography (abiotic conditions PC1 and PC2, and location). To minimize the influence of multicollinearity within each group (especially between abiotic covariates and location in the biogeography group; Supplementary Fig. 2a<\/a>), we performed variance partitioning at the level of variable groups. This was done by summing the variance terms of all predictors within each group81<\/a>. This allowed us to incorporate the covariances between group predictors to the variance component of that group, thereby reducing the covariance between variance components at the variable group level. Variance components of each group were standardized by the sum of all group-level variance components, including the residuals. To assess how covariance between group-level variance components might affect total variance accounting81<\/a>, we computed further variance partitioning metrics: (1) basic variance partition, which expresses the variance of each group standardized by the total variance of the dependent variable; (2) marginal variance partition, which quantifies the contribution of each group while accounting for covariances with all other groups; and (3) partial variance partition, which quantifies the contribution of each group that cannot be explained by a linear combination of the remaining groups. Further methodological details are provided in the Supplementary Methods<\/a>.<\/p>\n We calculated the effect size of fixed factors as the slope coefficient standardized by the standard deviation (Supplementary Methods<\/a>). We quantified the uncertainty of effect size based on credible intervals. Standardized slope (\u03b2st) values of |\u03b2st|\u2009<\u20090.12, 0.12\u2009\u2264\u2009|\u03b2st|\u2009<\u20090.24, 0.24\u2009\u2264\u2009|\u03b2st|\u2009<\u20090.41 and |\u03b2st|\u2009\u2265\u20090.41 were interpreted as neutral, weak, moderate and strong effects, respectively. We used P values of slope coefficients and Bayes factors as measures of evidence for fixed and random effects, respectively. The P values were derived from posterior distributions using a two-tailed test82<\/a>. Bayes factors were used to quantify the support for models including the random factor tested, compared with null models without the random factor. Bayes factors were calculated as the ratio of marginal likelihoods of the two models (Supplementary Methods<\/a>).<\/p>\n To test a priori multivariate causal hypotheses regarding how tree community effects on soil food web multifunctionality are mediated by changes in tree and understorey plant properties and microenvironment83<\/a>,84<\/a>, we performed piecewise multi-level structural equation modelling\u00a0(SEM)85<\/a>. We first constructed a causal diagram based on ecological theory and our previous knowledge of the system85<\/a> (Extended Data Fig. 1c<\/a> and Supplementary Fig. 3<\/a>) and built an initial SEM model (Supplementary Table 8<\/a>). Following a local estimation approach85<\/a>, we tested hypothesized causal relationships by fitting component LMMs with random intercepts across locations using Bayesian modelling, along with the use of variance partitioning and effect size computation methods in a manner similar to those described above. For each component LMM, abiotic covariates (see above) were included to ensure that mediation effects were not biased by confounding factors84<\/a>, thereby satisfying the \u2018backdoor criterion\u2019. The initial SEM model was then simplified by removing unsupported linkages85<\/a> using information-theoretic methods to select the most parsimonious SEM model. For each component LMM within the initial SEM model, stepwise backward selection was performed by removing the predictor with the highest P value, verifying at each step that model simplification improved goodness-of-fit using Pareto smoothed importance-sampling leave-one-out cross-validation86<\/a> (Supplementary Methods<\/a>). To test the goodness-of-fit of the selected SEM model, we used Shipley\u2019s test of directional separation87<\/a>, which assesses whether any important causal pathways may be missing (Supplementary Methods<\/a>). To quantify and compare the magnitude of tree community effects on soil food web multifunctionality, indirect effects were calculated as the product of the standardized coefficients along each path. Total effects were calculated as the sum of indirect effects.<\/p>\n Reporting summary<\/p>\n