Spatial patterns of soil multifunctionality trade-offs and their relationships with microbial diversity in warm-temperate mountain forests of North China

Soil ecological function (SEF) trade-offs are central to terrestrial ecosystem sustainability, yet their spatial patterns and relationships with microbial diversity in natural warm-temperate mountain forests remain poorly understood. Montane elevational gradients provide valuable natural settings for examining these relationships because marked hydrothermal variation over short distances, together with topsoil–subsoil differentiation, can reorganize microbial communities and the biogeochemical processes supporting multiple soil functions. We investigated three core SEFs—water and soil conservation (WSC), soil carbon sequestration stability (SCSS), and soil nutrient maintenance (SNM)—and their pairwise trade-offs along a 700–1,700 m elevational gradient in the topsoil and subsoil of Dongling Mountain, North China. High-throughput sequencing, multivariate statistical analyses, variation partitioning analysis (VPA), and Shapley additive explanations (SHAP) were used to evaluate the coupled responses of microbial diversity, individual SEFs, and SEF trade-offs. SEFs exhibited pronounced vertical stratification and elevational zonation. WSC and SNM were generally greater in the topsoil, whereas SCSS was higher in the subsoil. Among the three pairwise relationships, the SCSS–SNM trade-off was the strongest, reaching its maximum intensity at 1,300 m in the topsoil and 1,500 m in the subsoil. Bacterial and fungal diversity showed contrasting relationships with individual SEFs and nonlinear responses to trade-off intensity. VPA showed that altitude and soil depth explained the largest independent fraction of variation in individual SEFs, whereas soil properties explained the largest independent fraction of variation in SEF trade-offs. Among the microbial variables, bacterial diversity had the greatest relative predictive importance for individual SEFs and their trade-offs, while fungal diversity provided complementary predictive information. SHAP analysis further identified bacterial diversity as the most informative microbial predictor of overall soil multifunctionality and the soil N/P ratio as the strongest predictor of mean trade-off intensity. These findings indicate that soil multifunctionality and its trade-offs covary with soil depth, hydrothermal conditions, microbial diversity, and stoichiometric constraints. Within the sampled Dongling Mountain transect, this study provides quantitative evidence for spatially structured soil multifunctionality trade-offs and establishes a basis for testing the generality of these relationships across replicated mountain systems under global change.

Soil is a fundamental and nonrenewable component of terrestrial ecosystems and serves as a nexus of biogeochemical cycles, plant productivity, and environmental regulation (Philippot et al., 2023). Its multifunctionality, such as hydrological regulation (Li Y. et al., 2025), carbon sequestration (Wang S. Q. et al., 2025), nutrient cycling (Lu et al., 2025), and biodiversity maintenance (Wan et al., 2025), underpins global ecological balance, ecosystem stability, and resilience to global change (Lehmann et al., 2020). Previous studies of soil multifunctionality have largely focused on identifying its environmental drivers (e.g., climate, soil physicochemical properties, and vegetation) (Lehmann et al., 2020; Luo et al., 2025; Wang H. X. et al., 2025) and assessing its responses to anthropogenic and environmental disturbances (Liu H. Y. et al., 2025; Wang X. et al., 2025). However, these studies have predominantly treated soil multifunctionality as an integrated index or focused on individual soil ecological functions (SEFs) (He Y. et al., 2024; Xiao et al., 2025), largely overlooking the trade-offs and synergies among multiple SEFs that collectively determine the overall functioning of soil ecosystems.

Trade-offs among SEFs arise from inherent resource allocation constraints in natural soil systems, where the enhancement of one function often occurs at the expense of another (Luo et al., 2024). For example, the accumulation of stable organic matter to enhance carbon sequestration may reduce organic matter decomposition rates (Liang et al., 2025), thereby limiting nutrient mineralization and availability (Mei et al., 2024; Nan et al., 2025). Such trade-offs are ubiquitous in terrestrial ecosystems and directly modulate the long-term sustainability of soil ecosystem services (Daryanto et al., 2019). Montane elevational transects are particularly informative for examining these trade-offs because they encompass pronounced changes in hydrothermal and edaphic conditions over relatively short geographic distances (Körner, 2007; Sundqvist et al., 2013). In Dongling Mountain, the 700–1700 m transect spans mean annual temperatures of 4.2–8.1 °C and mean annual precipitation of 530–728 mm, with precipitation showing an overall increase with elevation (He L. B. et al., 2024a). Previous studies in this warm-temperate forest system further documented depth-dependent elevational structuring of bacterial communities and elevation- and depth-related shifts in fungal community composition and metabolic characteristics (He L.B. et al., 2024a,b). These coupled hydrothermal, edaphic, and microbial gradients are directly relevant to soil water retention, organic-matter transformation and stabilization, and nutrient cycling, which underpin water and soil conservation, soil carbon sequestration stability, and soil nutrient maintenance. Sampling both topsoil and subsoil adds a complementary vertical environmental axis, allowing us to examine whether the three SEFs, and the balance among them, co-vary across contrasting montane soil environments. We therefore treat elevation as an integrative observational gradient rather than as an isolated causal factor. Despite their ecological significance, studies on SEF trade-offs remain fragmented and inadequate. Existing research has focused primarily on grasslands, planted forests, or agricultural ecosystems (Han C. et al., 2025; Liang et al., 2025; Liu Y. J. et al., 2025), with limited attention given to mountain forest ecosystems—especially natural warm-temperate mountain forests characterized by complex topography, distinct elevational hydrothermal gradients, and vertical soil stratification. In these forests, how SEF trade-offs (e.g., between soil water conservation and carbon sequestration or nutrient maintenance and carbon stability) vary along the elevational gradient and soil depth gradients and what factors drive these variations remain poorly understood.

Soil microbial communities (bacteria and fungi) are core biological drivers of SEFs and mediate key processes such as organic matter decomposition (Husain et al., 2024), nutrient mineralization (Bai et al., 2020; Cao et al., 2024), and soil structure formation (Ye et al., 2025). Extensive research has documented links between microbial diversity and soil multifunctionality (Han et al., 2022; Zhou X. et al., 2024). For example, bacterial alpha diversity is positively correlated with soil nutrient cycling and carbon sequestration in Tibetan Plateau grasslands, with its response to elevational gradients directly mediating changes in soil multifunctionality (Chen et al., 2022); a study along European latitudinal gradients identified fungal diversity as a key driver of soil multifunctionality (Han X. et al., 2025). However, existing research has focused almost exclusively on microbial regulation of overall soil multifunctionality or individual SEFs (Wang C. et al., 2023), with a critical lack of investigation into how microbial diversity responds to and regulates trade-offs among multiple SEFs. Several key questions remain unresolved: Do bacterial and fungal diversity exhibit distinct response patterns to different pairwise SEF trade-offs? What are the mechanisms by which microbial diversity regulates the intensity of SEF trade-offs? These ambiguities stem from the complexity of microbial-SEF interactions (Wang et al., 2021) and the historical neglect of research on “SEF trade-offs” in favor of studies focusing on individual functions or overall multifunctionality (Wubs, 2024). While bacterial communities are often associated with rapid nutrient cycling and fungal communities with stable carbon sequestration (Naz et al., 2025; Qian et al., 2022), it remains unclear whether bacterial diversity decreases and fungal diversity increases (or vice versa, or weak correlations) when there is a trade-off between nutrient mineralization and carbon sequestration (Song et al., 2017).

Using this coupled elevational–vertical framework, we systematically investigated three core SEFs, their inter-function trade-offs, and their associations with soil microbial communities in the warm-temperate mountain forests of Dongling Mountain, North China. Soil samples were collected along a 700–1,700 m elevational gradient and from two soil depths, with three specific research objectives: (1) to characterize the spatial distribution patterns of the three core SEFs across the elevational gradient and vertical soil depths; (2) to quantify the intensity of trade-offs among SEFs and identify their environmental dependencies using the root mean square deviation (RMSD) method; and (3) to characterize the associations between soil bacterial and fungal diversity and the intensity of SEF trade-offs. We proposed two testable hypotheses: (1) SEFs and their trade-offs show significant elevational zonation and vertical stratification, with the strongest trade-off intensities occurring at mid-elevations; (2) bacterial and fungal diversity display divergent linear associations with individual SEFs and distinct nonlinear associations with SEF trade-off intensities, with bacterial diversity showing greater relative predictive importance than fungal diversity. Accordingly, this study provides a site-specific empirical basis for understanding SEF patterns and trade-offs along the Dongling Mountain transect and the development of sustainable ecosystem management strategies for warm-temperate mountain forests.

Study area and soil sample collection

This study was conducted in the Dongling Mountain region, which is located within Beijing Municipality, China. The area extends geographically from 39°48′ to 40°03′N and 115°24′ to 115°35′E, with the main peak reaching an elevation of 2,303 m. An elevational transect spanning 700–1,700 m above sea level was established along the mountain side. The region experiences a typical warm-temperate continental monsoon climate, with mean annual temperatures ranging from 4.2 °C to 8.1 °C and mean annual precipitation varying between 530 mm and 728 mm, with an increasing trend with elevation (He L. B. et al., 2024a). The zonal vegetation is a warm-temperate deciduous broadleaf forest, characterized by dominant tree species such as Populus davidiana, Betula platyphylla, and Quercus mongolica. To minimize major vegetation-related confounding, all 18 plots were confined to the contiguous 700–1,700 m warm-temperate deciduous broadleaf forest belt; areas below approximately 700 m, which have experienced strong anthropogenic disturbance including agricultural conversion, and the alpine meadow zone above 1,700 m were excluded. Vegetation type and dominant canopy taxa were recorded before soil sampling. Thus, the plots were consistent at the broad vegetation-formation level, although fine-scale variation in species abundance and stand structure could still occur along the transect (Supplementary Table S1).

Soil sampling was carried out in August 2022. Before the investigation and sampling, it was determined that there was no large-scale heavy rainfall for 7 days, which avoided short-term extreme climate interference. Six elevation gradients were selected: 700, 900, 1,100, 1,300, 1,500, and 1,700 m. The geographic location of the study area and the spatial arrangement of the six elevation-level sampling localities are shown in Supplementary Figure S1. At each elevation, three replicate plots of 20 m × 20 m were established within this forest formation on southeast-facing slopes with gradients <10° to reduce major variation in vegetation type, slope, and aspect among sampling locations. Within each plot, five soil cores were randomly collected from two depth intervals, 0–20 cm (topsoil) and 20–40 cm (subsoil), and combined into a single composite sample per depth. All samples were immediately transported to the laboratory on ice, sieved through a 2 mm mesh to remove visible plant roots and gravel, and then subdivided: one portion was stored at −80 °C for microbial analysis, another at 4 °C for measuring available nutrients, and the remaining portion was air-dried for assessing soil carbon, nitrogen content, and texture. Meteorological data, including mean annual temperature and precipitation for each sampling site, were obtained from the Resource and Environment Science and Data Center of the Chinese Academy of Sciences1 and the National Earth System Science Data Center.2

Selection and determination of soil ecological function indicators

On the basis of the collected soil sample data, we quantified 12 soil indicators spanning three categories of soil ecological functions (SEFs). Water and soil conservation: soil water storage (SWS) and soil erodibility index (SEI); soil carbon sequestration stability: dissolved organic carbon (DOC), free particulate organic carbon (fPOC), occluded particulate organic carbon (oPOC), mineral-associated organic carbon (MAOC), and microbial carbon (MBC); and soil nutrient maintenance: total nitrogen (TN), total phosphorus (TP), available phosphorus (AP), total potassium (TK), and available potassium (AK). Detailed information is provided in Supplementary Table S2.

The entropy weight method was employed to determine the weights of SEF indicators and conduct a comprehensive assessment of SEFs (You et al., 2025). The normalization, entropy, weighting, and aggregation procedures are given in Equations 1–6:

B ij and R ij denote the standardized and original measurements, respectively; R max and R min represent the maximum and minimum values, respectively, of the same indicator across different evaluation objects; P ij represents the contribution degree of the jth evaluation indicator to the ith evaluation object; E j refers to the entropy value of the jth evaluation indicator; W j denotes the weight of each evaluation indicator; and n represents the total number of indicators. SEFD is the comprehensive score for soil ecological function indicators.

Assessment of trade-offs in soil ecological functions

We quantified the trade-offs between SEFs using the root mean square deviation (RMSD) method (Ding et al., 2025), using Equation 7: where RMSD represents the root mean square deviation between two SEFs and is interpreted here as an index of relative imbalance among standardized functional values; E i is the normalized value for SEF i, calculated Equations 1 and 2; and E exp is the mean value across m SEFs.

To facilitate interpretation, each pair of normalized SEF values was projected in a two-dimensional coordinate system (Figure 1), following the geometric interpretation of RMSD used in ecosystem-service trade-off analyses (Lu et al., 2014; Ding et al., 2025). The 1:1 line (y = x) represents equal standardized values of the two functions, whereas the perpendicular departure of a paired observation from this line reflects the RMSD-based magnitude of imbalance; the side of the line indicates which SEF has the larger relative value. According to the RMSD results, we divided the trade-offs into five levels: high synergy (RMSD < 0.10), mild synergy (0.10 ≤ RMSD < 0.20), mild trade-off (0.20 ≤ RMSD < 0.30), medium trade-off (0.30 ≤ RMSD < 0.40), and high trade-off (0.40 ≤ RMSD). These categories are used only as operational descriptors of RMSD-based functional imbalance; in particular, the terms “high synergy” and “mild synergy” denote low imbalance and should not be interpreted as independent statistical evidence of positive covariation or causal facilitation. A detailed explanation of the schematic elements is provided in the Figure 1 legend.

DNA extraction, PCR amplification, high-throughput sequencing, and bioinformatics

Genomic DNA was extracted from 0.5 g of fresh soil samples (stored at −80 °C) using the MP FastDNA Spin Kit for Soil (Omega Biotek, Norcross, GA, U.S.) following the manufacturer’s standard protocols. The quality and concentration of the extracted DNA were verified by 1% agarose gel electrophoresis and a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, U.S.), ensuring that the A260/A280 ratio ranged between 1.8 and 2.0 for subsequent experiments. For bacterial community analysis, the V3-V4 hypervariable region of the 16S rRNA gene was amplified using the universal primer pair 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). For fungal community analysis, the ITS1 region of the rRNA gene was targeted with the primer pair ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2 (5′-GCTGCGTTCTTCATCGATGC-3′). PCR amplification was performed in a 25 μL reaction system containing 12.5 μL of 2 × Taq Plus Master Mix (Vazyme, Nanjing, China), 1 μL of each forward and reverse primer (10 μM), 2 μL of genomic DNA template (10 ng/μL), and 8.5 μL of sterile double-distilled water. The thermal cycling program for PCR was as follows: initial denaturation at 95 °C for 3 min; 35 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C (for bacteria) or 52 °C (for fungi) for 30 s, and extension at 72 °C for 45 s; final extension at 72 °C for 10 min, followed by storage at 4 °C.

PCR products were visualized by 1.5% agarose gel electrophoresis to check for target band specificity. Amplicons with clear and single bands were pooled in equal molar ratios and purified using the GeneJET Gel Extraction Kit (Thermo Scientific, Waltham, MA, U.S.). The purified products were quantified using a Qubit 3.0 fluorometer (Invitrogen, Carlsbad, CA, U.S.) and then subjected to high-throughput sequencing on the Illumina MiSeq PE300 platform (Illumina Inc., San Diego, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). The soil bacterial and fungal diversity data used in this study were derived from the 0–20 cm and 20–40 cm layers at elevations of 700–1,700 m on Dongling Mountain. These microbial sequencing data have been partially reported in two previously published studies by our research team (He L. B. et al., 2024a; He L. B. et al., 2024b), with the original sequencing data deposited in the NCBI Sequence Read Archive (SRA) under accession numbers PRJNA976495 (bacteria) and PRJNA1110162 (fungi).

Paired-end bacterial and fungal reads were processed using the dataset-specific workflows described in Supplementary methods S1 and clustered into operational taxonomic units (OTUs) at 97% sequence similarity. After chimera removal, bacterial and fungal OTUs were taxonomically assigned against the SILVA 16S rRNA and UNITE ITS reference databases, respectively. Shannon diversity, Chao1 richness, and Pielou evenness were calculated at the OTU level, whereas Bray–Curtis dissimilarities were used for NMDS and ANOSIM, with ANOSIM significance assessed using 999 permutations for both the bacterial and fungal datasets. Detailed quality-control procedures, software and database versions, sample-specific sequencing depths, sequencing-depth handling, and diversity calculations are provided in Supplementary methods S1 and Supplementary Table S7.

Data analyses and statistics

Before parametric analyses, the normality of the ANOVA residuals was assessed using the Shapiro–Wilk test, and homogeneity of variances was evaluated using Levene’s test. Box–Cox sensitivity analyses for variables showing departures from normality did not alter the significance patterns. Two-way ANOVA was used to test variation associated with elevation class, soil depth, and their interaction within the Dongling Mountain transect. Differences among elevations within each soil layer were compared using Duncan’s multiple comparison test at p < 0.05. Pearson correlation analysis was used to assess the relationships between microbial diversity indices and SEF indicators. Redundancy analysis (RDA) was performed to characterize associations of environmental variables and microbial diversity with SEFs and SEF trade-offs. Variance partitioning analysis (VPA) was further applied to quantify the independent and shared contributions of altitude-depth, soil properties, and microbial diversity to variations in SEFs and trade-offs. Linear regression was used to assess the relationships between microbial diversity and individual SEFs, whereas quadratic polynomial regression was used to characterize nonlinear relationships between microbial diversity and SEF trade-off intensity. Random forest (RF) and SHAP analyses were used to rank and interpret predictors of SEFs, soil multifunctionality, and trade-offs; RF importance was expressed as %IncMSE, and model tuning used shuffled five-fold cross-validation (seed = 42). Two-way ANOVA, Duncan’s multiple comparison test, Pearson correlation analysis, and linear and quadratic polynomial regression were performed in OriginPro 2026 (OriginLab Corporation, Northampton, MA, USA). Redundancy analysis and variance partitioning analysis were conducted in Canoco 5.0 (Microcomputer Power, Ithaca, NY, USA). The RF and SHAP analyses were implemented in Python 3.10 using the scikit-learn (v1.2.0) and shap (v0.41.0) packages. Given the modest sample size, these outputs were treated as exploratory within the sampled environmental range.

Distribution characteristics of SEF indicators

SEF indicators exhibited significant vertical and elevational heterogeneity (Figure 2). Across the 700–1,700 m gradient, SWS increased by 78.82% in the topsoil and 88.17% in the subsoil, although the topsoil response was nonmonotonic; both layers reached their maximum SWS at 1,700 m. The SEI was significantly lowest at 900 m in both layers (p < 0.05) and then generally increased with elevation. DOC was significantly greater at 900–1,300 m than at 1,500–1,700 m in both soil layers (p < 0.05); averaged across elevations, subsoil DOC was 30.49% lower than topsoil DOC. In the topsoil, fPOC and oPOC increased by 45.82 and 42.36%, respectively, from 700 to 1,700 m, whereas the subsoil fractions fluctuated without a corresponding monotonic increase and subsoil oPOC decreased by 80.87% between the two elevation of 700 m and 1,700 m. Across the entire gradient, mean subsoil fPOC and oPOC were 69.02 and 76.58% lower, respectively, than their topsoil values. MAOC decreased initially and then increased in both layers; from the 900 m minimum to the 1,700 m maximum, it increased by 86.85% in the topsoil, while those in the subsoil increased by a factor of 1.97. MBC showed contrasting altitude responses, decreasing by 26.82% in the topsoil but increasing by 43.28% in the subsoil from 700 to 1,700 m; nevertheless, mean subsoil MBC remained 76.37% lower than mean topsoil MBC across the gradient. TN in the topsoil remained comparatively stable and was highest at 1,700 m, whereas subsoil TN varied markedly and peaked at 1,300 m (p < 0.05). TP was highest in the topsoil at 1,700 m (p < 0.05) and in the subsoil at 1,300 m (p < 0.05). AP reached its maximum in the topsoil at 1,500 m and in the subsoil at 1,300 m, whereas TK in both layers increased to 1,100 m and then decreased. AK fluctuated but showed contrasting overall directions, decreasing by 28.55% between 700 and 1,700 m in the topsoil and increasing by 61.96% in the subsoil. Two-way ANOVA showed that elevation significantly affected SWS, SEI, DOC, MAOC, MBC, TN, TP, and AP, whereas soil depth significantly affected SEI, DOC, fPOC, oPOC, MAOC, MBC, TN, TP, TK, and AK (p < 0.05; Figure 2 and Supplementary Table S3). The elevation × depth interaction was significant only for oPOC (F = 4.620, p = 0.004) and MBC (F = 8.901, p < 0.001). Thus, depth effects were widespread, underscoring the importance of soil depth as a key factor driving variability in the SEF indicators.

Distribution characteristics of SEFs

On the basis of the 12 SEF indicators, we employed the entropy weighting method to determine the relative weights for the three SEFs and the 12 indicators, thereby establishing a comprehensive evaluation system for warm-temperate mountain forest soils (Supplementary Table S4). The three SEFs and soil multifunctionality (SMF) differed significantly across elevations and soil layers (Figure 3). WSC exhibited a clear hump-shaped elevational pattern (Figure 3a), with consistently higher values in the topsoil than in the subsoil. In both soil layers, WSC peaked at 900 m and was significantly higher at this elevation than at the other elevations (p < 0.05). In the topsoil, WSC increased from 0.04 at 700 m to 0.11 at 900 m, before declining to 0.05 at 1,700 m. The subsoil followed a similar pattern but showed a smaller amplitude. SCSS exhibited the opposite vertical pattern (Figure 3b), with mean values across elevations being 37.27% higher in the subsoil than in the topsoil. Topsoil SCSS showed greater elevational variation, increasing by 71.05% from its minimum at 1,300 m to its maximum at 1,700 m. By contrast, subsoil SCSS varied within a narrower range and peaked at 1,500 m, while the pooled 0–40 cm soil profile remained relatively stable across the elevational gradient. SNM was consistently greater in the topsoil (Figure 3c), with its across-elevation mean 91.79% higher than that of the subsoil; both layers displayed a mid-elevation maximum at 1,300 m followed by a decline at 1,500 m. SMF was slightly lower in the topsoil only at 1,500 m and higher at the other elevations (Figure 3d), yielding an across-elevation topsoil mean 15.41% higher than the subsoil mean. In the subsoil, SMF increased from 0.32 at 700 m to 0.42 at 1,300 m and then declined to 0.37 at 1,700 m, demonstrating a distinct mid-elevation optimum. Moreover, compared with the other elevations, both the subsoil and the whole soil layer presented significantly greater SMF at 1,300 m (p < 0.05). Two-way ANOVA showed significant elevation and soil-depth effects on WSC, SCSS, SNM, and SMF (Supplementary Table S5, p < 0.05). SCSS showed the strongest responses to elevation (F = 23.362, p < 0.001) and depth (F = 196.652, p < 0.001), and it was the only function with a significant elevation × depth interaction (F = 5.945, p = 0.001).

Changes in SEF trade-offs

We quantified the trade-offs among SEFs using the root mean square deviation (RMSD) method through pairwise analysis of WSC, SCSS, and SNM. As presented in Supplementary Table S5 and Figure 4, RMSD patterns depended on both elevation and soil depth. Soil depth significantly affected only the WSC–SCSS trade-off (F = 4.408, p = 0.046), whereas the elevation × depth interaction significantly affected both WSC–SCSS (F = 5.455, p = 0.002) and SCSS–SNM (F = 6.193, p < 0.001). The mean RMSD values across the three function pairs in the topsoil were 0.19, 0.31, 0.17, 0.32, 0.12, and 0.18 at 700, 900, 1,100, 1,300, 1,500, and 1,700 m, respectively. Thus, the topsoil maximum at 1,300 m was 168.93% greater than the minimum at 1,500 m; 1,300 and 900 m indicated medium trade-offs, whereas the other elevations showed mild synergy. In the subsoil, the corresponding mean RMSD values were 0.31, 0.17, 0.26, 0.22, 0.34, and 0.26, respectively. The subsoil maximum at 1,500 m was 98.12% greater than the 900 m minimum; 1,500 and 700 m showed medium trade-offs, 900 m fell within mild synergy, and the others showed mild trade-offs.

As shown in Figure 4, the diagonal plots for the different SEF pairs further elucidated the elevation- and depth-dependent trade-off patterns. For WSC–SCSS, the scatter points were widely distributed relative to the y = x equilibrium line (Figure 4a), indicating substantial variability in trade-off strength. In the topsoil, RMSD initially increased and then decreased with elevation, rising from 0.06 at 700 m (high synergy) to 0.47 at 900 m (high trade-off), with a significant difference between the two elevations (p < 0.05; Figure 4c). In the subsoil, RMSD initially decreased and subsequently increased, reaching its minimum at 900 m (0.16, mild synergy) and maximum at 1,500 m (0.41, high trade-off; Figure 4d). At the remaining subsoil elevations, WSC–SCSS trade-offs ranged from mild to high. In contrast, WSC–SNM points clustered near the y = x equilibrium line, indicating the weakest overall trade-off, with across-elevation mean RMSD values of 0.15 (mild synergy) in the topsoil and 0.12 (mild synergy) in the subsoil. This pair peaked at 1300 m in both layers (0.20 and 0.23, mild trade-off, respectively) and declined to 0.06 and 0.05 at 1700 m (high synergy), confirming that WSC and SNM generally remained comparatively balanced. Most SCSS–SNM points were located farther from the y = x equilibrium line, indicating stronger trade-offs. Its mean RMSD values were 0.25 (mild trade-off) in the topsoil and 0.33 (medium trade-off) in the subsoil. Topsoil SCSS–SNM RMSD peaked at 1,300 m (0.46, high trade-off) and then declined to 0.07–0.11 at 1,500–1,700 m, whereas subsoil RMSD reached its minimum at 1,300 m (0.13, mild synergy) and maximum at 1,500 m (0.52, high trade-off); both contrasts were significant (p < 0.05; Figures 4c,d). This reversal highlights a pronounced shift in the balance between carbon sequestration stability and nutrient maintenance across soil layers. Collectively, these results demonstrate that the aggregate elevational response reflected distinct, depth-specific trajectories among function pairs rather than a uniform change in all trade-offs.

Characteristics of soil microbial community composition and diversity

Nonmetric multidimensional scaling (NMDS) ordination based on Bray–Curtis dissimilarity revealed distinct elevational zonation of both bacterial and fungal communities in the topsoil and subsoil (p < 0.01) and low stress values (0.115 for bacteria and 0.145 for fungi), indicating a reliable ordination result (Figures 5a,d). Samples from the same elevation clustered into closely connected, independent groups along the NMDS1 axis. In contrast, the topsoil (T) and subsoil (S) samples at the same elevation showed substantial overlap within their respective elevation clusters. Circular taxonomic-composition chord diagrams further illustrated the phylum-level distribution of microbial communities across elevations and soil layers (Figures 5b,c,e,f). The bacterial profiles were dominated by Acidobacteria, Proteobacteria, and Chloroflexi (Figures 5b,c), with their relative representation varying among elevation groups. The fungal profiles were dominated by Basidiomycota and Ascomycota, with smaller contributions from Mortierellomycota and unclassified fungi (Figures 5e,f). Comparison of the topsoil (Figures 5b,e) and subsoil (Figures 5c,f) chord diagrams showed that the dominant phyla were broadly consistent between soil layers, whereas the relative representation of less abundant taxa differed, indicating depth-related compositional heterogeneity.

Two-way ANOVA indicated distinct responses of microbial α-diversity to elevation and soil depth (Figure 6; Supplementary Table S8). For bacteria, Shannon responded to elevation (F = 2.895, p = 0.035) and depth (F = 6.630, p = 0.017), Chao1 responded to elevation (F = 4.554, p = 0.005) and depth (F = 5.446, p = 0.028), whereas Pielou responded to depth (F = 7.867, p = 0.010) but not elevation (F = 1.796, p = 0.193). For fungi, depth affected Shannon (F = 5.935, p = 0.023), and both elevation (F = 2.641, p = 0.049) and depth (F = 115.250, p < 0.001) affected Chao1; Pielou was unaffected by elevation or depth (p ≥ 0.428). Consistent with these results, compared with the subsoil, most microbial α-diversity metrics were higher in the topsoil across the 700–1,700 m gradient. The bacterial Shannon and Pielou indices peaked at low elevations (700–900 m) and decreased at 1500 m (p < 0.05). However, the fungal Chao1 index was consistently higher in the topsoil than in the subsoil at every elevation, with mean values of 837.18 and 472.65, respectively. In the subsoil, fungal Chao1 increased by 83.01% from its minimum at 900 m to its maximum at 1,700 m (p < 0.05). In contrast, the fungal Shannon and Pielou indices did not significantly increase. The vertical divergence in bacterial diversity was most pronounced at mid-elevations (900–1,300 m), whereas fungal Chao1 richness was consistently stratified by soil depth regardless of elevation.

Relationships among soil microbial diversity, environmental factors and SEFs, and SEF trade-offs

Pearson correlation analysis (Figure 7a; Supplementary Table S6) indicated that bacterial α-diversity indices were associated with a more restricted subset of SEF indicators than fungal indices, particularly fungal Chao1 richness. B-Shannon was positively correlated with DOC, MBC, and AK (p < 0.05), B-Chao1 with SWS, fPOC, and MAOC (p < 0.05), and B-Pielou only with DOC (p < 0.05). In contrast, F-Shannon was positively correlated with fPOC, oPOC, MBC, TN, TP, and AK (p < 0.05), whereas F-Chao1 was positively correlated with all five SCSS indicators (DOC, fPOC, oPOC, MAOC, and MBC) and with TN, TP, and AK (p < 0.01). Notably, negative correlations between microbial diversity and SEF indicators were rare. RDA ordination (Figure 7b) was used to quantify the explanatory power of microbial diversity and environmental factors for multiple SEFs, with the first axis explaining 50.70% (topsoil) and 61.40% (subsoil) of the total variation. The associations of these factors with SEFs differed across soil layers: topsoil SEFs were most strongly associated with pH, B-Chao1, Pre, and Te, with pH and B-Chao1 positively associated with WSC, pH, B-Pielou, and Te negatively correlated with SCSS, and B-Pielou positively correlated with SNM (p < 0.05); subsoil SEFs were most strongly associated with pH, C/P, C/N, and Te, with B-Chao1 and F-Chao1 negatively associated with WSC, C/N positively correlated with SCSS, and C/P, N/P positively correlated with SNM (p < 0.05). The RDA ordination of the SEF trade-offs (Figure 7c) revealed that the first axis explained 50.32% (topsoil) and 63.03% (subsoil) of the total variation in the trade-offs, with the associated factors showing clear vertical differentiation. Topsoil SEF trade-offs were most strongly associated with B-Pielou, C/P, and EC; C/P was negatively associated with WSC-SCSS trade-off intensity; EC positively correlated with WSC-SNM intensity; and B-Shannon, B-Pielou, Te, and pH were positively associated with SCSS-SNM intensity (p < 0.05). Subsoil SEF trade-offs were most strongly associated with C/N, pH, B-Pielou, and F-Pielou; the pH was negatively correlated with the WSC-SCSS intensity, the B-Shannon index was negatively correlated with the WSC-SNM intensity, and C/N was positively associated with the SCSS-SNM intensity (p < 0.05). Variance partitioning analysis (VPA) further quantified the relative contributions of altitude-depth, soil properties, and microbial diversity to the variation in SEFs and SEF trade-offs. For SEFs, altitude-depth accounted for the largest independent fraction of the variance (27.44%), followed by microbial diversity (2.50%) and soil properties (2.13%), while their shared effects were also non-negligible, particularly the three-way shared fraction (21.09%); the residual fraction remained 36.33%. In contrast, for SEF trade-offs, soil properties explained the largest independent fraction (11.96%), exceeding altitude-depth (6.56%) and microbial diversity (2.68%), whereas the residual fraction increased to 62.38%. These results indicate that SEFs were primarily structured by topographic-vertical gradients, whereas SEF trade-offs were more strongly constrained by soil physicochemical properties, with substantial unexplained variation suggesting the involvement of additional unmeasured drivers.

Linear and quadratic polynomial regression analyses further revealed the response relationships between soil microbial diversity and multiple SEFs and SEF trade-offs (Figure 8). In terms of the linear relationships between soil microbial diversity and multiple SEFs (Figures 8a–c), bacterial diversity was positively associated with WSC (F = 4.263, p = 0.047) and SNM (F = 4.150, p = 0.049), whereas it was negatively associated with SCSS (F = 8.947, p = 0.005). The linear responses of fungal diversity displayed functional specificity: it was negatively associated with SCSS (F = 7.292, p = 0.011) and positively associated with SNM (F = 18.719, p < 0.001) but was not linearly associated with WSC (F = 1.672, p = 0.204). These results are consistent with the conclusions from the Pearson correlation analysis and redundancy analysis (RDA), confirming the different responses of bacteria and fungi to multiple SEFs. Quadratic polynomial regression further revealed the correlation characteristics between soil microbial diversity and SEF trade-offs (Figures 8d–f). With respect to the WSC-SNM trade-off, bacterial and fungal diversity exhibited significant “U-shaped” quadratic relationships (F = 6.299, p = 0.005 and F = 5.419, p = 0.009, respectively), indicating that microbial diversity was lowest at moderate trade-off intensities and significantly greater at both low and high trade-off intensities. With respect to the SCSS-SNM trade-off, bacterial diversity showed a significant “inverted U-shaped” quadratic relationship (F = 5.498, p = 0.009), with bacterial diversity peaking at moderate trade-off intensity and gradually decreasing at both ends; however, fungal diversity had no significant quadratic association with this trade-off (F = 0.403, p = 0.671). With respect to the WSC-SCSS trade-off, fungal diversity presented a significant “inverted U-shaped” quadratic relationship (F = 5.314, p = 0.010); fungal diversity was greatest at moderate trade-off intensities, whereas bacterial diversity was not significantly associated with this trade-off in the quadratic model (F = 1.912, p = 0.164). In summary, individual SEFs were linearly associated with microbial diversity, whereas SEF trade-off intensity showed nonlinear quadratic associations with microbial diversity. Bacterial diversity showed its strongest nonlinear associations with the WSC-SNM and SCSS-SNM trade-offs, whereas fungal diversity showed its strongest nonlinear associations with the WSC-SNM and WSC-SCSS trade-offs. These findings indicate distinct patterns of association in relation to soil multifunctionality.

Random forest (RF) analysis identified distinct dominant predictors of individual soil ecological functions (SEFs) and their trade-offs, highlighting strong functional specificity in predictor importance (Figure 9). For individual SEFs, bacterial diversity (B-Diversity), soil depth, and soil stoichiometric ratios were the primary predictors. Specifically, WSC was most strongly predicted by B-Diversity, soil depth, and N/P ratio (p < 0.05), whereas fungal diversity (F-Diversity) and hydrothermal factors contributed minimally. The strongest predictors of SCSS were soil depth and B-Diversity, both showing significant explanatory power (p < 0.05), while elevation, EC, and hydrothermal variables had relatively weak effects. SNM was strongly associated with soil depth, C/P ratio, B-Diversity, F-Diversity, and N/P ratio, with soil depth and C/P ratio acting as the dominant abiotic regulators. Trade-offs showed similarly differentiated predictor profiles. The WSC–SCSS trade-off was most strongly predicted by N/P ratio, B-Diversity, and EC, whereas the WSC–SNM trade-off was most strongly predicted by EC and N/P ratio, with microbial effects being comparatively weaker. In contrast, the SCSS–SNM trade-off was most strongly predicted by B-Diversity and N/P ratio and showed the strongest microbial dependence among all trade-offs, reflecting the close coupling between carbon stabilization and nutrient cycling processes. SHAP analysis further quantified the relative importance and directional effects of predictors on overall soil multifunctionality (SMF) and average trade-off intensity. B-Diversity exhibited the highest mean absolute SHAP value for SMF, indicating that it was the most informative microbial predictor. High B-Diversity was associated with greater predicted SMF, whereas low diversity was associated with lower predicted multifunctionality. Soil N/P ratio, elevation, and C/P ratio were secondary but important predictors, with high N/P and C/P ratios generally promoting SMF, while elevation exerted bidirectional effects. Soil depth, temperature, and precipitation had moderate influences, with topsoil consistently supporting greater SMF than subsoil. In contrast, average trade-off intensity was primarily controlled by N/P ratio, followed by C/N ratio and B-Diversity. High N/P and C/N ratios amplified trade-off intensity, whereas moderate B-Diversity corresponded to the lowest trade-off intensity, consistent with the nonlinear relationships identified in regression analyses. F-Diversity, pH, and C/P ratio also contributed substantially, with high fungal diversity and low pH generally intensifying SEF trade-offs.

Spatial distribution characteristics and environmental associations of soil ecological functions

The three SEFs showed clear vertical stratification and elevational zonation. WSC and SNM were consistently higher in the topsoil, whereas SCSS was higher in the subsoil. This pattern is likely related to the vertical distribution of soil organic matter and microbial activity. Topsoil receives continuous litter and root exudate inputs, which can stimulate microbial metabolism, nutrient cycling, and aggregate formation, thereby supporting water retention and nutrient maintenance (Huang W. et al., 2023; Zhou et al., 2022). In contrast, subsoil contains less labile organic matter but stronger mineral protection effects, which may favor the accumulation of stable carbon fractions such as MAOC and oPOC and, in turn, support SCSS (Li J. et al., 2023; Zhou Z. et al., 2024). Along the Dongling Mountain transect, elevation was associated with coordinated variation in hydrothermal conditions and potentially in vegetation and other soil-forming factors, which together may have shaped the observed SEF patterns. Within this integrated gradient, WSC peaked at mid-elevation, likely because increased precipitation enhanced soil water storage, while changes in vegetation composition and soil erodibility may have weakened soil structure at higher elevations (Guo et al., 2021; Wei et al., 2024). SCSS showed a U-shaped elevational trend, with the lowest values at 900 m, which may reflect the combined influence of temperature on organic matter decomposition and litter quality (Prescott and Vesterdal, 2021). SNM peaked at 1300 m, consistent with the commonly observed mid-elevation nutrient peak in mountain ecosystems (Huang C. et al., 2023; Zhang et al., 2023). Moderate hydrothermal conditions at mid-elevations may promote microbial mineralization while reducing nutrient loss through leaching (Bing et al., 2015). Two-way ANOVA and RDA further indicated that elevation, soil depth, stoichiometric ratios, and physicochemical properties were jointly associated with SEF variation. Soil stoichiometric ratios, particularly C/P and N/P, appeared to be important in linking microbial nutrient acquisition with carbon stabilization processes (Cui et al., 2022; Schleuss et al., 2021; Tan et al., 2025; Wu G. et al., 2023), suggesting that SEF heterogeneity along this transect was associated with interactions among hydrothermal gradients, vertical soil differentiation, and microbially mediated biogeochemical processes. Together, these depth-dependent responses across a spatially compressed environmental gradient support the use of a montane transect to examine variation in SEFs; however, this observational design requires treating elevation as a comprehensive environmental gradient rather than an isolated causal driver.

Cross-ecosystem comparisons indicate that the vertical differentiation observed in Dongling Mountain reflects a broadly applicable resource-gradient mechanism, whereas the precise elevational optima remain context dependent. Across soil profiles along temperate-forest elevation–bioclimate gradients, microbial community composition and extracellular-enzyme activity changed markedly with depth, while the direct effects of climate weakened and soil chemistry became relatively more important in deeper horizons (Dove et al., 2021). Depth-dependent responses of soil functionality and microbial communities have also been reported in subtropical forests, where tree-diversification effects differed among soil layers (Ding et al., 2024). Similarly, Alpine forest studies documented pronounced but site- and season-dependent shifts in soil nutrients, microbial abundance, and activity along elevation gradients (Siles et al., 2016), and work on the Tibetan Plateau showed that elevation-dependent microbial network complexity was closely associated with ecosystem functioning (Chen et al., 2022). Together, these studies support vertical resource differentiation—declining labile carbon and nutrient inputs, stronger mineral protection, and increasingly dominant edaphic filtering with depth—as a general mechanism. However, they do not support a universal elevation at which multifunctionality is maximized. Thus, the 900–1,300 m optima of WSC and SNM and the contrasting SCSS pattern observed here should be interpreted as transect-specific outcomes of locally covarying hydrothermal, vegetation, edaphic, and stoichiometric conditions rather than as general elevational thresholds.

Characteristics and formation mechanisms of soil ecological function trade-offs

Among the three pairs of trade-offs, the SCSS-SNM trade-off was the strongest, followed by the WSC-SCSS trade-off, while the WSC-SNM trade-off was the weakest (Figure 4). This difference in trade-off intensity is closely related to the resource demand and ecological processes of functional components in mountain forest ecosystems (Dantas de Paula et al., 2021). Soil carbon storage relies mainly on the accumulation of stable organic carbon fractions (e.g., MAOC and oPOC), which require the immobilization of nutrients such as nitrogen and phosphorus in soil organic matter (Cotrufo et al., 2019; Wu J. et al., 2023). Such nutrient immobilization may reduce nutrient bioavailability for plant uptake and microbial mineralization (Mao X. et al., 2024), providing a plausible basis for the comparatively high SCSS–SNM RMSD observed here. The WSC–SCSS trade-off may primarily depend on soil texture and structure. Soils with high clay contents can increase water retention capacity and simultaneously promote organic carbon stabilization through physical protection (Mao H.-R. et al., 2024); however, highly sandy soils may weaken both water retention capacity and carbon stabilization (Huang and Hartemink, 2020), resulting in a moderate trade-off between these two ecological functions. The weak WSC-SNM trade-off may stem from the synergistic promotion of these two ecological functions by the input of forest soil organic matter. Forest litter layers and root exudates not only increase the stability of soil aggregates and water-holding capacity but also provide a continuous nutrient supply for soil nutrient cycling (Chen et al., 2025; Panchal et al., 2022), thus offsetting their potential trade-off effects.

SEF trade-off patterns also varied with soil depth and elevation class along the Dongling Mountain transect. In the topsoil, the SCSS-SNM trade-off was the most significant at the mid-elevation of 1,300 m (Figure 4c), which is attributed to the high organic matter content and vigorous microbial activity in this elevational zone. Microbial decomposition of labile organic matter releases nutrients, which promotes soil nutrient cycling and storage (Huang and Su, 2025; Siegwart et al., 2023) while weakening the accumulation of stable organic carbon (Wu G. P. et al., 2024), thus exacerbating the trade-off between them. In the subsoil, the SCSS-SNM trade-off was the most significant at the higher elevation of 1,500 m (Figure 4d), which may be related to the insufficient nutrient supply in high-altitude areas (Cao et al., 2025). Soil microorganisms compete with plants for limited nutrients, leading to the prioritization of soil nutrient mineralization over organic carbon sequestration (Čapek et al., 2021; Wu H. et al., 2024), thus strengthening the SCSS-SNM trade-off. The WSC–SCSS trade-off differed significantly in the vertical direction: in the topsoil, this paired trade-off was the most significant at mid-low elevation (900 m) (Figure 4c), whereas in the subsoil, it was the most significant in the high-altitude zone (1,100–1,700 m) (Figure 4d), which was associated mainly by the distribution pattern of soil organic carbon fractions. In the topsoil, the contents of labile organic carbon such as DOC and fPOC are high at mid–low elevations (Figure 2). These substances may migrate downward with soil runoff at high altitudes, increasing soil water content while weakening soil carbon sequestration through leaching (Liang et al., 2021; Wang Z. et al., 2023); in the subsoil, the accumulation of stable organic carbon at high altitudes improves carbon sequestration capacity (Wu et al., 2025) but has little effect on water storage, thus resulting in a significant WSC–SCSS trade-off. These results highlight the strong environmental dependence of SEF trade-off patterns within the sampled Dongling Mountain transect.

Placing the SEF trade-off patterns within broader biodiversity–ecosystem-functioning, resource-allocation, and ecological-stoichiometry frameworks further clarifies their theoretical significance. Across global drylands and contrasting Scottish land-use systems, microbial diversity was positively related to an aggregate multifunctionality index even after accounting for climate, soil, and spatial covariates (Delgado-Baquerizo et al., 2016). Experimental grassland microcosms similarly showed that bacterial and fungal richness, functional redundancy, and interkingdom network complexity jointly supported nutrient-cycling functions (Wagg et al., 2019). These studies show that diverse microbial communities can maintain several functions simultaneously, but aggregate indices can conceal unequal or opposing responses among component functions. Our pairwise RMSD analysis adds this missing balance dimension. The comparatively small WSC-SNM imbalance is consistent with shared benefits from litter and root inputs, whereas the stronger SCSS-SNM imbalance is consistent with a finite-resource constraint: microbial processing and nutrient immobilization can promote oPOC/MAOC formation, while nutrient mineralization can increase carbon turnover, with the direction and magnitude conditioned by C:N:P stoichiometry (Cotrufo et al., 2019; Schleuss et al., 2021). By explicitly separating oPOC/MAOC stabilization from nutrient maintenance, the present study therefore extends work based only on bulk soil carbon or composite ecosystem-service indices. The ordering SCSS-SNM > WSC-SCSS > WSC-SNM is specific to the functions and environmental range sampled here, whereas the underlying principle - context-dependent allocation of finite water, carbon substrates, and nutrients among simultaneously demanded functions - is transferable across ecosystems.

Differentiated links between microbial diversity, soil ecological functions, and their trade-offs

Soil bacterial and fungal diversity exhibited different response patterns to SEFs and their RMSD-based relationships, reflecting distinct associations with multiple soil ecological processes. Both bacterial and fungal diversity were negatively correlated with SCSS and positively correlated with SNM, whereas only bacterial diversity was positively correlated with WSC (Figures 7, 8). These patterns may be related to the contrasting ecological roles of bacteria and fungi. Bacteria are commonly involved in rapid nutrient cycling and the degradation of labile organic matter (Raza et al., 2023). Greater bacterial diversity may therefore be associated with nutrient mineralization and with the formation of soil microaggregates through microbial biomass and extracellular polymeric substances (Costa et al., 2018; Gong et al., 2024; Yang Y. et al., 2023), processes that are relevant to SNM and WSC. The negative association between bacterial diversity and SCSS may reflect greater turnover of labile organic carbon under communities with broader bacterial resource-use capacities (Huang et al., 2021), although the present data do not identify the direction of causality. Fungi are often associated with the decomposition of recalcitrant organic matter and macroaggregate formation (Niego et al., 2023). Their positive association with SNM may reflect the gradual release of nutrients during decomposition of structurally complex substrates (Xue et al., 2023), whereas the negative association with SCSS may indicate that fungal-mediated decomposition in these warm-temperate soils did not translate directly into accumulation of highly stable mineral-associated carbon (Mayer et al., 2021). The absence of a significant linear relationship between fungal diversity and WSC suggests that, within this transect, water-holding capacity was more closely associated with bacterial diversity and abiotic soil properties than with fungal α-diversity alone.

Notably, microbial diversity showed distinct nonlinear responses to SEF trade-off intensities, revealing the adaptive strategies of microbial communities to functional trade-offs (Figure 8). Both bacterial and fungal diversity exhibited a significant “U-shaped” quadratic relationship with the WSC-SNM trade-off (p < 0.05), indicating that microbial diversity was the lowest at moderate trade-off intensities and significantly greater at both low and high trade-off intensities. At a moderate trade-off intensity between WSC and SNM, the soil ecosystem is in a relatively balanced state, the resource allocation between water and soil conservation and nutrient maintenance is optimized, and the microbial community is not subject to strong selection pressure, thus maintaining a moderate level of diversity; at a low trade-off intensity (synergy), sufficient organic matter input and superior hydrothermal conditions drive the synergistic enhancement of WSC and SNM, thus supporting high microbial diversity (Li F. et al., 2023; Zhang et al., 2025); at a high trade-off intensity, intense resource competition between WSC and SNM results in strong environmental selection pressure, promoting the survival of highly adaptable specialist microbial taxa (Zhang et al., 2024) and thus increasing microbial diversity. Bacterial diversity showed a significant “inverted U-shaped” quadratic relationship with the SCSS-SNM trade-off (p < 0.05), with bacterial diversity peaking at moderate trade-off intensities and gradually decreasing at both ends. A moderate SCSS-SNM trade-off intensity implies moderate resource competition between soil carbon sequestration and nutrient maintenance, a state that supports the coexistence of carbon sequestration-specialized and nutrient cycling-specialized bacterial groups (Hao et al., 2026; Piton et al., 2023), thus achieving the highest bacterial diversity; at a low trade-off intensity (synergy), the dominant process is the synchronous enhancement of SCSS and SNM, which results in the selection of a narrow range of microbial taxa adapted to this synergistic environment (Bao et al., 2023; Yang X. et al., 2023), thus reducing bacterial diversity; at a high trade-off intensity, intense resource competition leads to the dominance of a few highly competitive microbial taxa (Morrison et al., 2022; Wang and Kuzyakov, 2024), also reducing bacterial diversity. Fungal diversity exhibited a significant “inverted U-shaped” quadratic relationship with the WSC-SCSS trade-off (p < 0.05), indicating that fungal diversity peaked at moderate trade-off intensities. This pattern occurs because a moderate WSC-SCSS trade-off intensity can balance the demands of water retention and carbon stabilization on fungal community composition, allowing the coexistence of fungi adapted to high soil moisture and those specialized in carbon sequestration (Basiru and Hijri, 2024; Niu et al., 2024); at low or high trade-off intensities, strong environmental filtering leads to the dominance of fungi from a single functional group, resulting in a decline in fungal diversity (Glassman et al., 2017). These findings complement previous studies focusing on linear correlations between microbial diversity and SMF, and reveal the nonlinear response mechanism of microorganisms to SEF trade-offs.

Random forest analysis identified bacterial diversity as the strongest microbial predictor of SEFs and their trade-offs, whereas fungal diversity played a complementary role (Figure 9). Bacterial diversity showed stronger predictive power for WSC, SCSS, and most trade-offs, while fungal diversity showed stronger predictive associations with SCSS and the SCSS-SNM trade-off. This functional differentiation likely reflects the ecological roles of bacteria and fungi. Bacteria are key regulators of rapid nutrient cycling and soil structure formation (Devi et al., 2024), thereby influencing multiple SEFs and trade-offs, whereas fungi primarily contribute to the decomposition of recalcitrant organic matter and long-term carbon stabilization (Moravcová et al., 2025; Ullah et al., 2023). In addition, bacterial communities were more sensitive to elevational and vertical environmental gradients, whereas fungal communities exhibited greater environmental adaptability. Microbial diversity also interacted strongly with abiotic factors, including soil stoichiometric ratios, pH, and electrical conductivity, forming a coupled “abiotic–microbial–function” association structure (Figures 7, 9). For example, the SCSS-SNM trade-off was jointly associated with bacterial diversity and soil N/P ratio, whereas the strongest predictors of the WSC-SNM trade-off were bacterial diversity, N/P ratio, and EC. These findings highlight the necessity of simultaneously considering biotic and abiotic drivers when evaluating SEF trade-offs. SHAP analysis further quantified the relative importance and directional effects of predictors (Figure 9). Bacterial diversity exhibited the highest mean absolute SHAP value for overall SMF, indicating that it was the most important microbial predictor in this site’s environment. High bacterial diversity was associated with higher predicted SMF, whereas low diversity was associated with lower predicted SMF. Soil N/P ratio was the strongest predictor of average trade-off intensity, with high N/ p values associated with higher predicted trade-off intensity. Notably, bacterial diversity showed a bidirectional predictive pattern for trade-off intensity, with moderate diversity corresponding to the lowest trade-off intensity, consistent with the nonlinear relationships identified by regression analysis.

Cross-ecosystem evidence also shows that the relative importance of bacterial and fungal diversity is context dependent. Our result is consistent with a subtropical forest study in which bacterial rather than fungal diversity and community assembly were more strongly associated with soil multifunctionality (Han et al., 2022), and with an alpine-grassland restoration experiment in which bacterial diversity and network stability were the strongest microbial predictors of multifunctionality (Shu et al., 2023). However, boreal-forest soils showed the opposite pattern, with fungal richness and community composition being more informative than bacterial richness (Li et al., 2019). Moreover, drought–rewetting responses in temperate forest soils differed between bacterial and fungal communities and depended on historical land use (Osburn et al., 2021), illustrating that disturbance history can alter which microbial group tracks ecosystem functioning. These supporting and contrasting findings accord with context-dependent biodiversity–ecosystem-functioning theory: complementarity and functional redundancy may strengthen multifunctionality, but environmental filtering, substrate quality, and disturbance can shift the relative contribution of bacterial and fungal assemblages. Therefore, bacterial α-diversity should be interpreted as the strongest microbial predictor for the sampled WSC–SCSS–SNM suite and environmental range of the Dongling Mountain transect, not as a universal or experimentally demonstrated regulator.

Ecological implications and future research directions

Based on a single transect at Dongling Mountain, this study offers site-specific management insights, while still providing a valuable theoretical and practical foundation for sustainable soil management in warm-temperate mountain forests. First, the spatial patterns of SEFs and their trade-offs indicate that management priorities should be soil-layer specific, because the SCSS-SNM trade-off peaked at 1,300 m in the topsoil and 1,500 m in the subsoil (Figure 4). Accordingly, protecting native forest vegetation, reducing human disturbances (e.g., logging and grazing), and maintaining stable organic carbon pools are relevant candidate measures for balancing SCSS and SNM. Second, the functional differentiation of microbial diversity suggests that maintaining bacterial and fungal diversity may be relevant to soil multifunctionality. Management measures such as maintaining organic inputs, avoiding excessive soil disturbance, and regulating soil pH may therefore be considered candidate directions, but the present observational analyses do not demonstrate that increasing microbial diversity will directly enhance specific SEFs. These associations identify candidate management directions that require validation before application beyond the study transect. Third, the observed associations with soil depth, elevation class, and soil stoichiometric ratios provide site-specific hypotheses about how SEF trade-off patterns may respond to global change. Future climate warming may alter the vertical distribution pattern of hydrothermal conditions, thereby changing the intensity of SEF trade-offs and the regulatory role of microbial diversity, and targeted adaptive management measures can be formulated accordingly.

Several limitations should be acknowledged. First, this study focused mainly on microbial α diversity, while functional genes, microbial interactions, and community assembly processes were not considered. Second, only three core SEFs were included in the trade-off analysis. Third, all six elevation levels were sampled along a single Dongling Mountain transect. The three plots at each elevation are local replicates within an elevation band and do not provide independent replication of elevation across separate mountains or transects. Although this design avoided major transitions in vegetation formation and land use, fine-scale variation in species abundance, stand structure, litter inputs, and root distribution may still have covaried with elevation. Because these vegetation attributes were not included as independent quantitative predictors, their contributions cannot be separated from the integrated elevation. Fourth, the regression, random forest, and SHAP analyses quantify associations and predictive importance, not causal regulation. Fifth, the RMSD-based trade-off measure used here quantifies relative imbalance among standardized SEF values. Although it is useful for comparing functional balance across elevations and soil depths, RMSD alone does not establish the direction of statistical covariation or demonstrate causal competition or facilitation between functions. Future studies should incorporate additional functions, microbial functional traits, replicated mountain systems, and manipulative or longitudinal designs. Within these constraints, the present study documents depth- and elevation-dependent variation in WSC, SCSS, SNM, and their relative imbalances, and identifies bacterial α-diversity as the strongest microbial predictor for this functional suite within the sampled environmental range. These findings provide new insights into microbial associations with soil multifunctionality trade-offs and provide a site-specific empirical basis for future replicated mountain studies and management evaluation.

This study systematically characterized the spatial patterns of three core SEFs and their pairwise trade-offs, together with their relationships to bacterial and fungal diversity along the 700–1,700 m Dongling Mountain elevational transect in North China. Our results demonstrate pronounced vertical stratification and elevational zonation of SEFs: WSC and SNM were generally greater in the 0–20 cm topsoil, whereas SCSS was higher in the 20–40 cm subsoil. Among the three paired functions, the SCSS–SNM trade-off was the strongest, with peak intensities at 1,300 m in the topsoil and 1,500 m in the subsoil. Microbial communities also varied systematically with elevation and soil depth, with consistently higher bacterial and fungal α diversity in the topsoil than in the subsoil across the full elevational gradient. Importantly, microbial diversity was not only associated with SEFs but was also linked to the intensity and direction of their trade-offs. Bacterial and fungal diversity showed different linear relationships with individual SEFs and distinct nonlinear responses to SEF trade-off intensity. Bacterial diversity was the most informative microbial predictor of overall soil multifunctionality, whereas fungal diversity appeared to play a complementary role. SHAP analysis further suggested that soil N/P ratio was the primary predictor of average SEF trade-off intensity, and moderate bacterial diversity was associated with the lowest trade-off intensity. These biotic patterns were jointly linked with abiotic factors, including soil depth, elevational hydrothermal conditions, and soil stoichiometric ratios. These transect-level findings fill a site-specific gap in SEF trade-off research and provide a basis for future replicated studies and locally evaluated management strategies under global change.