The study area is located in Xishuangbanna, Yunnan Province, China (21°10′–22°40′ N and 99°55′–101°50′ E), belonging to the tropical humid zone south of the Tropic of Cancer (Zhai et al., 2017). The climate in this region is characterized as tropical monsoon climate, with an average annual temperature ranging from 15.1 to 21.7°C and an annual rainfall between 1200 and 2500 millimeters (mm) (Cao et al., 2006). The region is affected by the southwest monsoon throughout the year, with pronounced dry and wet seasons. The soil type is categorized as acidic brick-red soil (Ling et al., 2022). This region has the largest remaining area of tropical rainforest in China. In our study region, the main tree species in the tropical rainforest included Parashorea chinensis, Cleidion brevipetiolatum, Phoebe lanceolata, and Garcinia cowa. The shrub species included Boehmeria glomerulifera and Pittosporopsis kerrii, while the liana species included Tetrastigma planicaule and Aganosma breviloba.

From April to May 2021, we established three 1 hm² forest dynamic plots in the Menglun, Shangyong, and Bubeng areas, totaling 75 sampling plots with the area of 20 m × 20 m. In each sampling plot, we identified and recorded all woody plants with a diameter at breast height (DBH) ≥ 1 cm and their tree height (m), including tree, liana, and shrub. The volume of plant branches after removing the bark was determined using the drainage method, and then the branches were dried in an oven at 103°C to a constant weight that was measured after drying. The branch density was calculated using the ratio of branch dry weight to volume. We calculated the above-ground biomass (AGB) among the woody plant species using species-specific allometric equations (Equation 1) based on wood density, DBH, and tree height (Byrnes et al., 2014; Li et al., 2020c).

where AGB represents the above-ground biomass of the plant, WD represents wood density, DBH represents diameter at breast height, and H represents tree height.

We collected a total of 75 soil samples. In each sample plot, we took five soil samples of about 250g from the center and four corners, with a 0–20 cm soil depth, and mixed them to form a composite soil sample. Then, we divided the composite soil sample into two parts for measuring the physicochemical properties of the soil and the soil microbial DNA, respectively. At the same sampling point, a stainless-steel cylinder with a volume of 100 cm³ was used to collect soil core, which were used to measure physical indicators such as soil bulk density (SBD), soil water content (SWC) and soil water-holding capacity (SWHC) (Reynolds, 1970; Wang et al., 2022b). Next, we used standard measurement methods to determine soil chemical indicators such as soil organic carbon (SOC), soil total nitrogen (TN), soil total phosphorus (TP), soil hydrolyzed nitrogen (HN), soil available phosphorus (AP), soil total potassium (TK), and soil available potassium (AK). We used the CHN analyzer (2400 II CHN Analyzer, PerkinElmer, Boston, MA, USA) to measure SOC and TN, and the molybdenum antimonial blue colorimetry and Olsen method to determine TP and AP (Bao, 2000; Li et al., 2022a). In addition, HN was determined by alkaline hydrolysis diffusion method. TK and AK were determined by atomic absorption spectrophotometry (Bao, 2000). The formula for calculating soil organic carbon stock (SOCS) (Equation 2) was used as follows:

where SOCS (t•hm^-2) represents soil organic carbon stock, c (g•kg^-1) represents soil organic carbon content, d (g•cm^-3) represents soil bulk density, r (cm) represents soil layer thickness, and D represents the gravel content with a diameter greater than 2 mm, measured as a volume ratio.

We extracted DNA from 0.5 g soil samples using the FastDNA Soil DNA Isolation Kit (MP Biochemicals, Solon, OH, USA). Soil microbial diversity was analyzed using the nucleotide conserved region of rRNA. For bacteria, we used the V4 - V5 region of the highly conserved 16S gene, and for fungi, we mainly used the internal transcriptional spacer (ITS) region. The bacterial 16S rRNA was amplified with gene primers 338F (5’-ACTCCTACGGGAGGCAGCA-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’). The ITS of fungi was amplified by primers ITS1 (5’-CTTGGTCATTTAGAGGAAGTAA-3’) and ITS2 (5’-GCTGCGTTCTTCATCGATGC-3’). After PCR, the amplicons were extracted and purified for downstream library preparation for high-throughput sequencing. We used the Illumina HiSeq sequencing platform to analyze soil microbial diversity and constructed short fragment libraries using paired-end sequencing methods to analyze microbial species annotation and abundance. Among them, polymerase chain reaction (PCR) was performed in triplicate using 50 μl of reaction mixture containing 10 μl buffer, 0.2 μl Q5 High-Fidelity DNA polymerase, 10 μl High GC Enhancer, 1 μl Dntp, 10 μM per primer, and 60 ng genomic DNA. The PCR for fungi was performed using the Veriti Thermal cycler (Applied Biosystems). The thermal cycle conditions were 95°C for 5 min, then 95°C for 1 min, 50°C for 1 min, 72°C for 1 min, a total of 15 cycles, and finally extended for 10 min at 72°C (Li et al., 2020c). We performed amplicon sequencing using the Illumina MixSeq 2500 platform of Biomarker Technologies Corporation, Beijing, China. Low quality sequences (quality score< 20, length< 150 bp, total expected error > 0.5) were filtered after combining forward and backward reads. After sequencing, we clustered or de-noised high-quality sequences of soil bacteria and fungi and divided the bacterial and fungal communities into operational taxonomic units (OTUs) based on 97% similarity. After sequencing 75 samples of soil bacteria and fungi, we obtained 6,001,239 and 6,001,299 pairs of end-reads, respectively, from which 79,167 and 79,221 clean reads were generated. The original sequencing data can be found in the NCBI Sequence Read Archive (SRA) under the accession number SUB13691156: PRJNA995948.

In this study, EMF include carbon stock, nutrient cycling, water regulation and biomass production. Among them, biomass production represented the above-ground biomass (AGB) of all woody plants with DBH ≥ 1 cm. At the same time, nutrient cycling was composed of soil total nitrogen (TN), soil total phosphorus (TP), soil hydrolyzed nitrogen (HN), soil available phosphorus (AP), soil total potassium (TK), and soil available potassium (AK) in the ecosystem. Soil organic carbon stocks (SOCS) served as a proxy for carbon stock, and soil water-holding capacity (SWHC) represented water regulation.

We standardized these four ecosystem functions across 75 sampling plots using Z-score transformations. Then, the four ecosystem function scores for each sampling plots were averaged to obtain a multifunctional index (Byrnes et al., 2014; Yang et al., 2023).

All of the statistical analyses were performed using R4.2.3 (R Core Team, 2022). First, we use a multithreshold multifunctional method to evaluate the effects of plant diversity and soil microbial diversity on EMF, and analyze the trade-offs between different ecosystem functions of EMF to verify whether plant diversity and soil microbial diversity support multiple ecosystem functions at high threshold levels. First, we run this method using the “multifunc” package. The effects of plant diversity and soil microbial diversity on EMF were examined by the slope changes of curves fitted with the EMF index (threshold-based index of multifunctionality (MFt)) at different thresholds (Byrnes et al., 2014). Among them, the minimum threshold when diversity has an effect on multifunctionality (Tmin), the maximum threshold when diversity has an insignificant effect on multifunctionality (Tmax), the threshold when diversity has the greatest positive or negative effect (Tmde), and the threshold when diversity has the greatest effect on the number of functions that reach the threshold (Rmde). These four main indicators give us information on how plant diversity and soil microbial diversity affect EMF.

Next, we used the “agricolae” package to conduct ordinary least squares regression analysis to explore the relationship between plant diversity (including tree, shrub, and liana) and soil microbial diversity, soil microbial network complexity, soil water content, soil bulk density, and EMF (R Core Team, 2022). Then, we used the “linkET” package to perform correlated analyses and visualizations of the relationships between plant diversity, soil microbial diversity, soil microbial network complexity, environmental factors, and EMF. We used the ‘WGCNA’ package to construct a co-occurrence network of soil bacteria and fungi based on Spearman correlation coefficient for microbial co-occurrence network analysis. In our study, only OTUs with a relative abundance > 0.01% were used in the microbial co-occurrence network, where OTUs represented nodes and the correlations between nodes represented edges (Qiu et al., 2021). We used the “RMThreshold” package to set appropriate similarity thresholds of 0.850 and 0.540 for soil bacteria and fungi based on random matrix theory (RMT). Then, we adjusted the P-values (P< 0.05) of the correlations using the Benjamini-Hochberg false discovery rate (FDR) (Benjamini et al., 2006). We use the “subgraph” function in the “igraph” package to extract the subnetwork properties of each soil sample from the constructed network (Yang et al., 2022). Then, the soil microbial network was visualized using the Gephi platform (Xue et al., 2017). Soil bacterial and fungal network complexity was calculated using the topological features of the network, including the number of nodes and edges, average degree, clustering coefficient, average path length, network diameter, and graph density (Jiao et al., 2022). We also checked the dominant taxonomic groups at the phylum level (Qiu et al., 2021).

Next, we use multidimensional scales (MDS) to calculate topological parameters of soil microbial network complexity. Before calculation, the average path length and diameter are processed using the reciprocal of the variable (X-1) (Jiao et al., 2022). Then, before statistical analysis, all data were log10-transformed for normality. OTUs represented the alpha diversity index of soil bacteria and fungi, and plant diversity is mainly represented by plant species richness, included the number of tree, shrub, and liana species in each sampling plot. The effects of plant diveristy, soil microbial diversity, soil microbial network complexity, and environmental factors on EMF were assessed using ordinary least squares linear regression.

Finally, we use “lavaan” (Rosseel, 2012) and “randomForest” packages to analyze SEM and random forest model respectively. We used random forest model mean square error (MSE) percentages to assess the relative importance of plant diversity, soil microbial diversity, soil microbial network complexity, and environmental factors to EMF (Jiao et al., 2022). Next, we used structural equation modeling (SEM) to assess the direct and indirect associations between plant diversity, soil microbial diversity, soil microbial network complexity, environmental factors, and EMF. The effects of different variables on EMF and multiple ecosystem functions were determined by path standardization coefficients and correlation P-values. We use Chi-square tests, goodness of fit index (GFI), and approximate root-mean-square error (RMSEA) to evaluate the fit of the model. The results showed that P > 0.05, GFI = 1, RMSEA< 0.05 by chi-square test (Rosseel, 2012), indicating that the structural equation model was trustworthy.

According to the results obtained from the multithreshold method, we can conclude that TSR has a significant positive correlation with EMF in the 1% to 3% threshold range ( Figure 1A ). LSR has a significant negative effect on EMF, with a minimum threshold of 13% and a maximum threshold of 81%. When the threshold reached 48%, the effect on EMF was the greatest, and the maximum effect (Rmde) of LSR was 0.1668, meaning that an increase of one liana species reduced 0.1668 functional scores ( Figure 1B ). SSR has no significant effect on EMF ( Figure 1C ). In addition, BSR has a significant positive effect on EMF, with a minimum threshold of 86% and a maximum threshold of 92% ( Figure 1D ). In the 3% to 60% threshold range, FSR has a significant negative effect on EMF ( Figure 1E ).

Based on ordinary least squares linear regression analysis, EMF decreased with increasing LSR ( Figure 2B ), BSR ( Figure 2D ), FSR ( Figure 2E ) and FNC ( Figure 2G ). However, EMF showed an increasing trend with the increase in BNC ( Figure 2F ). Additionally, there was no significant regression relationship between EMF and TSR ( Figure 2A ), SSR ( Figure 2C ), SWC ( Figure 2H ), and SBD ( Figure 2I ).

Further, combined heatmap analysis showed that FSR, FNC, and BNC had significant positive effects on nutrient cycling function, and the affected soil nutrients mainly included TN, TP, TK, HN, and AP ( Figure 3A ). SSR, TSR, and LSR mainly have significant positive effects on AP, and LSR also has significant positive effects on TN and SOCS ( Figure 3B ). SBD and SWC have significant positive effects on SWHC, and SWC also has a positive effect on AP ( Figure 3C ).

Soil bacterial and fungal communities in the soil exhibited differential patterns of symbiosis, but both soil bacterial and fungal co-occurrence networks demonstrated high complexity. Based on the topological parameters of the soil microbial networks ( Figure 4E ), the average degree of soil bacteria network was 16.889. The network diameter was 8, with a graph density of 0.158. The modularity was 0.599, and the average clustering coefficient of the nodes was 0.661, with an average path length of 2.647. On the other hand, the average degree of soil fungi network was 88.041. The network diameter was 5, with a graph density of 0.454. The modularity was 0.359, and the average clustering coefficient of the nodes was 0.786, with an average path length of 1.713.

At the phylum level, the dominant groups of soil bacteria and their proportion in the total phylum abundance ( Figures 4A, C ), were Acidobacteria (45.37%), Proteobacteria (24.07%), and Verrucomicrobia (11.11%), among others. The dominant groups of soil fungi were Ascomycota (52.57%) and Basidiomycota (35.43%) ( Figures 4B, D ). Among the dominant bacterial groups, Acidobacteria and Verrucomicrobia showed a significant negative correlation with EMF, mainly influencing the soil nutrient cycling function of TN, TP, TK, and SOCS. On the other hand, Proteobacteria showed a significant positive correlation with EMF, positively impacting the nutrient cycling of TN, TP, and TK and carbon stock function ( Figure 5 ). Among the dominant fungal groups, Ascomycota had a significant negative correlation with EMF, primarily affecting nutrient cycling function of AP, TP, and TK. Conversely, Basidiomycota were positively correlated with EMF, mainly influencing nutrient cycling function of TP, TK, and carbon stock function ( Figure 5 ).

According to the random forest model ( Figure 6B ), LSR, BNC, and FNC ranked the top three in the importance of contributing EMF, with importance scores of 15.37, 14.13, and 12.91, respectively, while FSR, SSR, and SWC ranked the least in the importance of EMF, with scores of 1.41, 0.74, and 0.64, respectively.

The standardized effects of plant diveristy, environmental factors, soil microbial diversity, and soil microbial network complexity on EMF were obtained through structural equation modeling (SEM) ( Figure 6A ). Chi-square test showed no significant P-value (0.926 > 0.05), indicating that there was no significant difference between the implied covariance matrix of the model and the observed covariance matrix (P< 0.05), which confirmed the validity of the structural equation model. The results demonstrated that BNC has a highly significant direct positive effect on EMF (path standardized coefficient = 0.471, P< 0.001). TSR also has a significant direct positive effect on EMF (path standardized coefficient = 0.306, P< 0.01), while LSR has a highly significant direct negative effect on EMF (path standardized coefficient = -0.450, P< 0.001). Furthermore, MSR primarily exerted an indirect effect on EMF by influencing MNC, while environmental factors (EF) had significant negative effects on LSR and FNC, thus indirectly affecting EMF.