Our study was conducted in the Wuyi Mountain National Nature Reserve (27°73′–27°86′N, 117°69′–117°78′E), located in the northwestern region of Fujian Province, China (Suppl. material 1: table S1). The reserve is characterized by a typical subtropical monsoon climate, marked by a pronounced altitudinal gradient and notable climatic variations. The mean annual temperature (MAT) ranges from 9.14°C to 14.65°C, and the mean annual precipitation (MAP) varies between 1961.00 mm and 2308.00 mm. Spanning an area of 290 square kilometers, the reserve is home to extensive native vegetation (Li et al. 2021) (Suppl. material 1: table S2). Five study sites were selected at elevations of 675 m, 1227 m, 1762 m, 1852 m, and 2157 m, representing distinct vegetation types: evergreen broad-leaved forest (EBF), coniferous and broad-leaved mixed forest (CBMF), coniferous forest (CF), subalpine dwarf forest (SDF), and alpine meadow (ALM) (Suppl. material 1: table S1). The dominant species in each respective forest type include Castanopsis eyrei, Syzygium austrosinense, and Cyclobalanopsis myrsinifolia in the EBF; Pinus taiwanensis, Rhododendron ovatum, Eurya muricata, and Cyclobalanopsis multinervis in the CBMF; Pinus taiwanensis and Tsuga chinensis in the CF; Eurya saxicola, Symplocos lucida, Yushania hirticaulis, Deyeuxia pyramidalis, and Cyperus microiria in the SDF; and Miscanthus sinensis, Aniselytron treutleri, and Carex olivacea in the ALM.

At each altitudinal level, three replicate study plots (20 m × 20 m) were systematically established. Each plot was subsequently divided into 16 contiguous quadrats (5 m × 5 m) through grid-based partitioning, with division boundaries precisely aligned at 5-meter intervals along both the × and y axes. A comprehensive tree survey was carried out within each quadrat. Three quadrats were randomly selected, and within each, five soil samples were collected from a depth of 0–15 cm in an S-shaped pattern. These samples were then homogenized and mixed to form a single composite sample. The samples were placed in sterile polyethylene bags, stored in an icebox, and transported to the laboratory within one day. Each composite sample was split into two subsamples: one for analysis of soil physicochemical properties and enzyme activities, and the other for storage at -80°C for total DNA extraction. The geographic coordinates of each sampling point were recorded using a global positioning system. Leaf Area Index (LAI), Enhanced Vegetation Index (EVI), and Normalized Difference Vegetation Index (NDVI) data were obtained from the LAADS DAAC (https://ladsweb.modaps.eosdis.nasa.gov) and extracted by latitude and longitude using the Environment for Visualizing Images software (Suppl. material 1: table S3).

The concentrations of soil total carbon (TC) and total nitrogen (TN) were analyzed using an elemental analyzer (FLASH SMART; Shanghai, China). Total phosphorus (TP) concentration was determined by alkali fusion molybdenum-antimony resistance spectrophotometry (Duan et al. 2023). Ammonium nitrogen (NH₄⁺-N) and nitrate nitrogen (NO₃⁻-N) concentrations were measured using indophenol blue colorimetry and phenol-disulfonic acid colorimetry, respectively (Zhao et al. 2022). Soil organic matter content was quantified through potassium dichromate oxidation spectrophotometry (Cen et al. 2021). Microbial biomass carbon (MBC), nitrogen (MBN), and phosphorus (MBP) were determined using the chloroform fumigation extraction method (Xu et al. 2021). Soil pH was measured with a pH meter (PB-10, Shanghai, China) in a 2.5:1 (water:soil) suspension, and electrical conductivity (EC) was measured with a conductivity meter (DDS-307A, Shanghai, China) in a 1:1.5 (water:soil) suspension. Soil moisture content was determined by the oven drying method (O’Kelly 2004).

The activities of soil enzymes were assessed as follows: urease activity was measured using indigo colorimetry, acid phosphatase activity by 2,3,5-triphenyl tetrazolium chloride colorimetry, dehydrogenase activity by triphenyltetrazolium chloride colorimetry, and leucine aminopeptidase activity by p-nitroaniline colorimetry. Cellulase and sucrase activities were determined using 3,5-dinitrosalicylic acid colorimetry, while β-1,4-glucosidase, β-1,4-N-acetylglucosaminidase, and β-xylosidase activities were quantified using the p-nitrophenol method (Mao and Jiang 2021; Xie et al. 2021; Boteva et al. 2022; Wang et al. 2022a; Lu et al. 2023; Qu et al. 2023).

Total DNA was extracted from soil samples using the Fast DNA™ Spin Kit for Soil (MP Biomedicals, California, USA) following the manufacturer’s protocol. The concentration and quality of the extracted DNA were evaluated by 1% agarose gel electrophoresis. For fungal community analysis, the ITS1 region of the rRNA gene was amplified using the primer pair ITS1F (CTTGGTCATTTAGAGGAAGTAA) and ITS2R (GCTGCGTTCTTCATCGATGC) (Bell et al. 2020). Polymerase chain reaction (PCR) was conducted in a 20-μl reaction mixture containing 10 μl of 2× Pro Taq Master Mix, 0.8 μl each of forward and reverse primers (5 μM), 50 ng of template DNA, and ddH₂O to adjust the final volume to 20 μl. The PCR cycling protocol included an initial denaturation at 95°C for 3 minutes, followed by 35 cycles of denaturation at 95°C for 30 seconds, annealing at 55°C for 30 seconds, and extension at 72°C for 45 seconds, with a final elongation step at 72°C for 10 minutes and a hold at 10°C.

For AMF, the V4–V5 hypervariable regions of the 18S rRNA gene were amplified using a nested PCR approach. The first round of PCR utilized the primer pair AML1F (ATCAACTTTCGATGGTAGGATAGA) and AML2R (GAACCCAAACACTTTGGTTTCC) (Yang et al. 2021) in a 20-μl reaction system with 2× Pro Taq Master Mix. The cycling conditions consisted of an initial denaturation at 95°C for 3 minutes, followed by 32 cycles of 95°C for 30 seconds, 55°C for 30 seconds, and 72°C for 45 seconds, with a final extension at 72°C for 10 minutes. The second round of PCR employed the primers AMV4-5NF (AAGCTCGTAGTTGAATTTCG) and AMDGR (CCCAACTATCCCTATTAATCAT) (Yang et al. 2021), with the following conditions: 95°C for 3 minutes, followed by 25 cycles of 95°C for 30 seconds, 55°C for 30 seconds, and 72°C for 45 seconds, and a final elongation at 72°C for 10 minutes. To prevent contamination, a negative control (ddH₂O) was included during DNA extraction and PCR procedures.

The PCR products were initially visualized on a 2% agarose gel and subsequently quantified using the QuantiFluor™-ST Blue fluorescence quantification system (Promega, Madison, USA). Based on the quantification results, the PCR products were pooled in equimolar ratios according to the sequencing requirements for each sample. Sequencing libraries were prepared using the TruSeq™ DNA Sample Prep Kit (Sangon Biotech, Shanghai, China) and sequenced on the Illumina HiSeq 2500 platform (Illumina, Inc., San Diego, USA) at Majorbio Co. (Shanghai, China).

Paired-end reads generated on the Illumina platform were first merged using FLASH version 1.2.11 (Magoc and Salzberg 2011) based on overlap relationships. Adapters and low-quality sequences were then trimmed using Trimmomatic version 0.33, and chimeric sequences were detected and removed with UCHIME version 8.1 (Edgar et al. 2011; Bolger et al. 2014). The resulting high-quality, non-chimeric sequences were clustered into operational taxonomic units (OTUs) at a 97% similarity threshold using USEARCH version 11 (Edgar 2013). Fungal taxonomic classification was performed by aligning sequences against the UNITE reference database (version 7.0; https://unite.ut.ee/) (Nilsson et al. 2019). AMF were taxonomically assigned by querying representative sequences against the MaarjAM database (Öpik et al. 2010), while EMF were identified by cross-referencing the FungalTraits database (version 1.2) for lineages and exploration types associated with ectomycorrhizal lifestyles (Põlme et al. 2020). OTU-level taxonomic assignments were conducted using the RDP Classifier version 2.13 with an 80% confidence threshold (Wang et al. 2007). To account for uneven sequencing depths across samples, sequences were rarefied to the minimum depth using the “sub.sample” command in MOTHUR (version 1.46.0) (Schloss et al. 2009).

Most data analyses in this study were performed using R (version 4.2.2). The α-diversity and phylogenetic diversity of mycorrhizal fungi were calculated using the “vegan” (version 2.7-0) and “picante” (version 1.8.2) R packages (Kembel et al. 2010). Fungal community differences were assessed through principal coordinates analysis (PCoA) and permutational multivariate analysis of variance (PERMANOVA). The influence of environmental factors on vertical community composition was evaluated using the “ecodist” (version 2.1.3) package (Goslee and Urban 2007). To explore the effects of altitude on fungal community composition, β-diversity was partitioned into species turnover and richness differences using the “adespatial” (version 0.3-24) package (Dray and Dufour 2007). Additionally, the “tidyverse” (version 2.0.0) package (Wickham et al. 2019) was employed to analyze the species composition of AMF and EMF in relation to altitudinal gradients.

Fungal community assembly mechanisms were elucidated by quantitatively evaluating the relative contributions of distinct ecological processes using the “iCAMP” package (version 1.5.12). In this framework, a βNRI value < -1.96 indicates homogeneous selection, while a βNRI value > 1.96 signifies heterogeneous selection. When |βNRI| ≤ 1.96 and the Raup-Crick (RC) value < -0.95, the process is identified as homogeneous dispersal. Conversely, |βNRI| ≤ 1.96 and RC > 0.95 suggest dispersal limitation. Finally, |βNRI| ≤ 1.96 and |RC| ≤ 0.95 are indicative of ecological drift (Ning et al. 2020). To identify potential keystone taxa within fungal communities, the specificity and occupancy rates of each OTU across altitudinal gradients were calculated and visualized (Gweon et al. 2021). The impact of environmental factors on the topological properties of mycorrhizal fungal networks was assessed using the “corr.test” function from the “psych” R package (version 2.4.6.26).

To investigate interactions within mycorrhizal fungal communities, co-occurrence networks were constructed using the weighted gene co-expression network analysis (WGCNA) framework in R (Langfelder and Horvath 2012). OTUs were represented as network nodes, and statistically significant associations between taxa were depicted as edges. To minimize false positives, multiple hypothesis testing corrections were applied using the “multtest” package (version 2.54.0) (Benjamini et al. 2006), with a significance threshold of α = 0.01. Indirect correlations were removed using network topology disentanglement algorithms (Feizi et al. 2015), retaining only robust pairwise relationships that met dual criteria: an absolute Spearman’s ρ > 0.6 and an adjusted P-value < 0.01 (Ma et al. 2015; Qian et al. 2019). Final network visualizations were produced using the Gephi platform (version 0.9.3; https://gephi.org), with node spatial arrangements optimized through force-directed algorithms (Mathieu et al. 2014).

Functional diversity metrics—including functional divergence and functional evenness—were calculated for mycorrhizal fungi based on community traits (e.g., richness, Shannon diversity, phylogenetic diversity, mean pairwise distance, network modularity, network density, node count, and edge count) and environmental preferences (e.g., pH, soil moisture content, MAT, and EC) using the “FD” package (version 1.0-12.3) (Escalas et al. 2019). Statistical significance across groups was assessed using one-way ANOVA and the aovMcomper function from the “EasyStat” package (version 0.1.0) when data met assumptions of normality and homogeneity of variance. Otherwise, the Kruskal-Wallis test and the KwWlx2 function were employed to evaluate overall and pairwise differences.

The influence of environmental factors on mycorrhizal fungal diversity was assessed using the Mantel test, implemented with the “linkET” package (version 0.0.7.3) (Liu et al. 2022). Following the approach of Yang et al. (2023b), TC, TN, TP, NO₃⁻-N, NH₄⁺-N, MBC, MBN, MBP, sucrase, urease, and leucine aminopeptidase were selected as proxies for ecological functions, representing soil nutrient storage, nutrient cycling, productivity, and fertility. Ecosystem multifunctionality was quantified using the mean value method, calculated as the average of Z-normalized individual function values. Structural equation modeling (SEM) and path coefficient analysis were conducted to evaluate the effects of altitude, climate, soil properties, and fungal communities on ecosystem multifunctionality, utilizing the “piecewiseSEM” package (version 2.3.0) (Liu et al. 2022). The relative importance of environmental factors to ecosystem multifunctionality was determined using a random forest model, implemented with the “rfPermute” package (version 2.5.2) (Zhang et al. 2022b). Non-linear regression was employed to model the relationships between mycorrhizal fungal diversity, altitude, and ecosystem multifunctionality, with the coefficient of determination (R²) and P-values of the fitted curves calculated using the lm function (Milici et al. 2016; Chen et al. 2022; Peng et al. 2023).

AMF were classified into five distinct orders, with Glomerales exhibiting the highest relative abundance, ranging from 79.14% to 97.70%. This order was most abundant in EBF at lower altitudes and least abundant in SDF at intermediate altitudes. Diversisporales followed, with a relative abundance ranging from 1.79% to 19.12% (Suppl. material 1: fig. S1A). In contrast, EMF were primarily divided into two phyla: Basidiomycota, which consistently exceeded 96% relative abundance across all altitudinal gradients, and Ascomycota, which showed a relative abundance ranging from 1.39% to 3.29% (Suppl. material 1: fig. S1B). At the genus level, the AMF community comprised 15 genera, while the EMF community included 62 genera. Glomus was the dominant genus for AMF across all altitudes. In contrast, the dominant EMF genera varied with altitude: Russula was prevalent in CBMF and SDF at intermediate altitudes; Sebacina dominated in EBF at the lowest altitudes; and Entoloma was most abundant in ALM at the highest altitudes (Suppl. material 1: fig. S1C, D).

AMF richness peaked significantly in SDF at intermediate altitudes (Fig. 1A). Within this community, Glomerales, the dominant order in Wuyi Mountain, exhibited no significant altitudinal variation in richness (Fig. 1B). In contrast, Diversisporales, the second most abundant order, showed considerable variability in richness across the studied altitudinal gradients (Fig. 1C). Similarly, the EMF community displayed significant variations in richness, with the highest values observed in EBF at the lowest altitude and the lowest in ALM at the highest altitude (Fig. 1D). Richness within the Basidiomycota and Ascomycota phyla of the ectomycorrhizal community also exhibited significant altitudinal variations (Fig. 1E, F). PCoA revealed distinct community structures for both AMF and EMF across different altitudes (Fig. 1G, H). This differentiation was further supported by PERMANOVA, which confirmed significant altitudinal effects on community composition for both fungal types (Fig. 1G, H; Suppl. material 1: table S4). Environmental factor analysis identified strong correlations between fungal community variations and a range of soil and climatic parameters, including soil moisture, enzymatic activities, microbial biomass, nutrient concentrations, and vegetation indices (Suppl. material 1: table S5). Decomposition of β-diversity indicated that species replacement was the primary driver of community differences in both AMF and EMF, accounting for 56.45% and 57.18% of the variation, respectively, with additional contributions from richness differences and similarity (Fig. 1I, J).

Figure 1. 

Altitudinal variation in mycorrhizal fungal richness and community differentiation. Richness of various taxa of mycorrhizal fungi including Arbuscular (A), Glomerales (B), Diversisporales (C), ectomycorrhizal fungi (D), Basidiomycota (E), and Ascomycota (F) across five altitudinal gradient zones. Community differences among arbuscular and ectomycorrhizal fungi are shown in G and H, respectively, while I and J depict the decomposition of β-diversity within these fungal communities. Repl: replacement; RichDiff: richness difference.

AMF richness showed no significant correlation with altitude (P = 0.601) (Suppl. material 1: fig. S2A). Glomerales exhibited consistent richness across altitudinal gradients (P = 0.218), whereas the richness of Diversisporales increased significantly with altitude (P < 0.001, R² = 0.481) (Suppl. material 1: fig. S2B, C). In contrast, EMF communities, including both Ascomycota and Basidiomycota, displayed significant altitudinal fluctuations in richness, characterized by an initial decline, followed by an increase, and a subsequent decrease (P < 0.05, R² > 0.230) (Suppl. material 1: fig. S2D–F). We further investigated altitudinal trends in phylogenetic and functional diversity among mycorrhizal fungi. AMF phylogenetic diversity and functional evenness showed no significant vertical variation (P > 0.560, R² = 0.048) (Suppl. material 1: fig. S2G, I). However, AMF functional divergence exhibited an inverted S-shaped pattern with increasing altitude (P < 0.001, R² = 0.566) (Suppl. material 1: fig. S2H). EMF phylogenetic diversity followed a similar fluctuating pattern with altitude, while functional evenness was inversely correlated, and functional divergence remained stable (P = 0.515, R² = 0.001) (Suppl. material 1: fig. S2G–L).

Further analyses examined the influence of altitude on mycorrhizal fungal community diversity and its associations with environmental variables. The richness of Diversisporales showed significant correlations with soil moisture content, MAT, precipitation, and the EVI (P < 0.05) (Suppl. material 1: fig. S3A). AMF functional divergence was significantly linked to TN, TP, organic matter, microbial biomass, soil moisture, climatic conditions, and sucrase activity (P < 0.05) (Fig. 2A). EMF richness, particularly within Basidiomycota, was significantly associated with TC, TN, TP, organic matter, microbial biomass, soil moisture, and climatic variables (P < 0.05) (Fig. 2B and Suppl. material 1: fig. S3B). Ascomycota richness was significantly related to MAT, precipitation, EVI, and dehydrogenase activity (P < 0.05) (Suppl. material 1: fig. S3C). EMF phylogenetic diversity was correlated with similar soil and climatic parameters, including urease activity (P < 0.05) (Fig. 2C). Finally, EMF functional evenness exhibited significant relationships with TP, MBN, MBP, soil moisture, climatic variables, urease, and leucine aminopeptidase (P < 0.05) (Fig. 2D).

Figure 2. 

Influence of environmental variables on mycorrhizal fungal diversity and function. Environmental factors impact the functional diversity of arbuscular mycorrhizal fungi (A), richness (B), phylogenetic diversity (C), and functional evenness (D) of ectomycorrhizal fungi. TC: total carbon; TN: total nitrogen; TP: total phosphorus; OM: organic matter; MBC: microbial biomass carbon; MBN: microbial biomass nitrogen; MBP: microbial biomass phosphorus; NO3._N: nitrate nitrogen; NH4._N: ammonium nitrogen; EC: electric conductivity; MC_F: fresh soil moisture content; MC_AD: air-dried soil moisture content; MAT: mean annual temperature; MAP: mean annual precipitation; CL: cellulase; UE: urease; SC: sucrase; ACP: acid phosphatase; β-GC: β-1,4-glucosidase; DHA: dehydrogenase; NAG: β-1,4-N-acetylglucosaminidase; β-XYS: β-xylosidase; LAP: Leucine aminopeptidase.

The iCAMP analysis revealed that dispersal limitation and ecological drift were the primary drivers structuring the assembly of the two distinct mycorrhizal fungal communities, collectively accounting for over 86.50% of the observed assembly processes (Suppl. material 1: fig. S4). In contrast, the contributions of the other three ecological processes were relatively minor, representing less than 13.45% of the assembly dynamics (Suppl. material 1: fig. S4). The assembly of mycorrhizal fungi was further influenced by both altitude and fungal type. Specifically, ecological drift dominated the community assembly of AMF, exhibiting a complex pattern of increase, decrease, and subsequent increase with altitude (Fig. 3A). In contrast, dispersal limitation was the primary factor shaping the assembly of EMF, displaying a pattern inverse to that of AMF (Fig. 3B). Homogenizing dispersal was observed exclusively at intermediate altitudes in AMF communities, while it occurred at both the lowest and intermediate altitudes in EMF communities (Fig. 3A, B). Heterogeneous selection was detected in both fungal types at the same altitude, though its influence varied between the communities (Fig. 3A, B).

Figure 3. 

Dominance of ecological processes in mycorrhizal fungal community assembly. The relative importance of heterogeneous selection, homogeneous selection, dispersal limitation, homogenizing dispersal, and drift in shaping arbuscular (A) and ectomycorrhizal (B) fungal communities. HeS: heterogeneous selection; HoS: homogeneous selection; DL: dispersal limitation; HD: homogenizing dispersal; DR: drift.

SPEC-OCCU plot analyses revealed that AMF and EMF OTUs across various altitudes did not cluster within specific x-axis columns, indicating a high diversity of occupancy characteristics at different elevations (Suppl. material 1: fig. S5). The study identified potential keystone taxa—defined as OTUs with specificity and occupancy values no less than 0.7—across five altitudinal stands, demarcated by dotted lines. As altitude increased, the distribution of potential keystone OTUs in AMF communities exhibited significant variation, representing 0.14, 0.33, 0.01, 0.27, and 0.28 of the community at each respective altitude (Suppl. material 1: fig. S5A–E). Similarly, the proportion of potential keystone OTUs within EMF communities varied, comprising 0.16, 0.65, 0.03, 0.56, and 0.62 from the lowest to the highest altitudes (Suppl. material 1: fig. S5F–J).

Within the AMF community, Glomerales were identified as potential keystone taxa across all forest belts. Diversisporales appeared at intermediate altitudes in CF and SDF, and at the highest altitudes in ALM, while Archaeosporales were exclusively observed in SDF at intermediate altitudes (Suppl. material 1: fig. S5A–E). For the EMF community, Basidiomycota and Ascomycota emerged as potential keystone taxa; Basidiomycota were found at all altitudes, while Ascomycota were specifically observed in the CBMF at intermediate altitudes and in ALM at the highest altitudes (Suppl. material 1: fig. S5F–J).

Co-occurrence networks for AMF and EMF consisted of 665 and 281 nodes, respectively, connected by 6,302 and 1,416 edges (Fig. 4A, G; Suppl. material 1: table S6). The AMF network exhibited greater complexity than the EMF network, although it displayed lower values for average clustering coefficient, average path length, diameter, density, and modularity (Fig. 4A, G; Suppl. material 1: table S6). Variations in network complexity were observed across altitudinal gradients. Specifically, the AMF network was most complex in EBF at the lowest altitude and least complex in ALM at the highest altitude. Conversely, the EMF network reached peak complexity in CF at intermediate altitudes and was least complex in ALM at the highest altitude (Fig. 4B–F, H–L; Suppl. material 1: table S6). Notably, no shared edges or nodes were observed between AMF and EMF networks across the five studied altitudes (Suppl. material 1: fig. S6). Topological features of these networks, including average clustering coefficient, path length, diameter, density, and modularity, exhibited significant altitudinal variations (Suppl. material 1: table S6).

Figure 4. 

Co-occurrence networks of mycorrhizal fungal communities across different ecosystems. Total and specific community networks of arbuscular and ectomycorrhizal fungi in various forest types (A–L), with analysis of the relationship between network topological features and environmental factors for arbuscular (M) and ectomycorrhizal fungi (N). EBF: evergreen broad-leaved forest; CBMF: coniferous and broad-leaved mixed forest; CF: coniferous forest; SDF: subalpine dwarf forest; ALM: alpine meadow. ACC: average clustering coefficient; APL: average path length.

Further analysis revealed relationships between network topology and environmental factors. In the AMF network, modularity was significantly associated with cellulase, urease, acid phosphatase, and soil moisture content. Network density showed strong correlations with urease, β-xylosidase, acid phosphatase, MBN, soil moisture content, and LAI. Diameter was significantly linked to dehydrogenase, MAT, precipitation, and EVI. Additionally, average path length correlated with TP, while nodes demonstrated significant associations with NO₃⁻-N and EC (Fig. 4M). In the EMF network, density was significantly correlated with NDVI, and the average clustering coefficient showed significant relationships with dehydrogenase, TP, soil moisture content, MAT, precipitation, and EVI. Significant associations were also observed between edges and NH₄⁺-N, EVI, while nodes were notably related to NDVI (Fig. 4N).

Ecosystem multifunctionality demonstrated significant variation along altitudinal gradients, with the highest levels observed in the ALM at the highest elevation and the lowest levels in the CBMF at intermediate elevations. To elucidate the drivers of this variability, we employed a SEM and analyzed the associated path coefficients. The results revealed both direct and indirect effects of multiple factors on multifunctionality. Soil properties emerged as the most influential factor, exhibiting a strong direct effect (0.697) and a moderate indirect effect (0.273). AMF also showed a substantial direct effect (0.697) (Fig. 5A, B). Climatic factors ranked third in terms of overall influence, with a total effect of 0.548 (Fig. 5A, B). Altitude had a direct positive effect (0.317) on multifunctionality but a slight negative indirect effect (-0.081), resulting in a net total effect of 0.236, which exceeded the effect of EMF alone (-0.127) (Fig. 5A, B). Additionally, random forest analysis (explained variation: 84.27%, P < 0.001) identified network modularity, network density, MAT, MAP, and altitude as the primary determinants of multifunctionality (Fig. 5C).

Figure 5. 

The relationships between ecosystem multifunctionality and environmental factors across altitudinal gradients. A, B direct and indirect effects of environmental factors on multifunctionality; C relative importance of individual environmental factors. EBF: evergreen broad-leaved forest; CBMF: coniferous and broad-leaved mixed forest; CF: coniferous forest; SDF: subalpine dwarf forest; ALM: alpine meadow; MC_F: fresh soil moisture content; MC_AD: air-dried soil moisture content; MAT: mean annual temperature; MAP: mean annual precipitation; AMF_PD; arbuscular mycorrhizal fungal phylogenetic diversity; AMF_rich: arbuscular mycorrhizal fungal richness; AMF_Modularity: arbuscular mycorrhizal fungal network modularity; AMF_Density: arbuscular mycorrhizal fungal network density; EMF_PD; ectomycorrhizal fungal phylogenetic diversity; EMF_rich: ectomycorrhizal fungal richness; EMF_Modularity: ectomycorrhizal fungal network modularity; EMF_Density: ectomycorrhizal fungal network density.

The authors have declared that no competing interests exist.

No ethical statement was reported.

All the fungal strains used in this study have been legally obtained, respecting the Convention on Biological Diversity (Rio Convention).

This work was supported by the Nature Science Foundation of Fujian province (2021J05023), the National Nature Science Foundation of China (32260027) and the Program for Young Talents of Science and Technology in Universities of Inner Mongolia Autonomous Region (NJYT24051).

Taotao Wei: Writing–original draft; Huiguang Zhang: Investigation; Shunfen Wang: Investigation, Formal analysis; Chunping Wu: Investigation; Tieyao Tu: Writing–review & editing, Investigation, Conceptualization; Yonglong Wang: Investigation, Methodology, Funding acquisition; Xin Qian: Writing–review & editing, Formal analysis, Investigation, Conceptualization, Methodology, Funding acquisition, Project administration.

Chunping Wu https://orcid.org/0009-0000-5888-5444

Data will be made available on request.