Linking Bacterial-Fungal Relationships to Microbial Diversity and Soil Nutrient Cycling

ABSTRACT Biodiversity is important for supporting ecosystem functioning. To evaluate the factors contributing to the strength of microbial diversity-function relationships in complex terrestrial ecosystems, we conducted a soil survey over different habitats, including an agricultural field, forest, wetland, grassland, and desert. Soil microbial multidiversity was estimated by the combination of bacterial and fungal diversity. Soil ecosystem functions were evaluated using a multinutrient cycling index (MNC) in relation to carbon, nitrate, phosphorus, and potassium cycling. Significant positive relationships between soil multidiversity and multinutrient cycling were observed in all habitats, except the grassland and desert. Specifically, community compositions showed stronger correlations with multinutrient cycling than α-diversity, indicating the crucial role of microbial community composition differences on soil nutrient cycling. Importantly, we revealed that changes in both the neutral processes (Sloan neutral modeling) and the proportion of negative bacterial-fungal associations were linked to the magnitude and direction of the diversity-MNC relationships. The habitats less governed by neutral processes and dominated by negative bacterial-fungal associations exhibited stronger negative microbial α-diversity–MNC relationships. Our findings suggested that the balance between positive and negative bacterial-fungal associations was connected to the link between soil biodiversity and ecosystem function in complex terrestrial ecosystems. This study elucidates the potential factors influencing diversity-function relationships, thereby enabling future studies to forecast the effects of belowground biodiversity on ecosystem function. IMPORTANCE The relationships between soil biodiversity and ecosystem functions are an important yet poorly understood topic in microbial ecology. This study presents an exploratory effort to gain predictive understanding of the factors driving the relationships between microbial diversity and potential soil nutrient cycling in complex terrestrial ecosystems. Our structural equation modeling and random forest analysis revealed that the balance between positive and negative bacterial-fungal associations was clearly linked to the strength of the relationships between soil microbial diversity and multiple nutrients cycling across different habitats. This study revealed the potential factors underpinning diversity-function relationships in terrestrial ecosystems and thus helps us to manage soil microbial communities for better provisioning of key ecosystem services.

relative importance of these processes changes in space and time (2)(3)(4). Deterministic processes involve nonrandom and niche-based mechanisms (5), including environmental filtering and interspecific interactions (e.g., competition, facilitation, mutualisms, and predation). In contrast, stochastic processes reflect mainly random changes in the relative abundances of species, involving random birth, death, and dispersal events (6,7). Uncovering the community assembly processes is critical to understanding the generation of biodiversity and its contribution to ecosystem functions (8)(9)(10). On one hand, stochastic assembly processes were dominant in high-a-diversity communities, associating with the existence of specialized functions that are correlated with specific bacterial taxa (9). On the other hand, the prediction of niche theory stating that competing species lead to reduced niche overlap implies a negative covariation between species through competition (11). For example, the global niche differentiation between fungi and bacteria relating to contrasting diversity responses to precipitation (12) might induce negative bacterial-fungal covariation, whereas niche partitioning may promote positive species covariations due to fitness differences among organisms under environmental heterogeneity and the niche-partitioning scenario (13,14), which might consequently increase functional community performance (15,16). For example, bacteria profiting from the organic matter degradation of fungi (17) may induce their positive covariation. Therefore, revealing the balance between deterministic and stochastic processes may help to predict the functional contributions of ecological communities, which should be higher in real-world environments than in simplified experimental settings (11). Given this point, it is of great importance to understand how these naturally diverse and fluctuating communities are organized (e.g., community assembly) and how they influence the functioning of ecosystems.
Biodiversity is known as a critical determinant of ecosystem functioning (18). Understanding the relationship between biodiversity and ecosystem functioning is a central issue in ecology that contributes to better provision of key ecosystem services to humans (19). The strength of the biodiversity-ecosystem function (BEF) relationship may be determined by functional redundancy, complementarity, or species competition (19)(20)(21). For example, losses of one or a few species in a community with high functional redundancy may have minimal consequences for ecosystem function (22); the complementarity effect describes that niche differentiation or facilitation between species can increase the functional performance of communities (23), and species competition for resources can reduce community functioning when species performing major functions are inhibited (20). However, these previous studies have been based on well-manipulated experimental conditions (16,19,20); most natural communities feature a highly complex taxonomic diversity (18,24). Currently, it is very challenging to characterize the determinants of BEF relationships in complex natural ecosystems.
The activities of microorganisms and their interactions greatly influence a variety of ecosystem processes associating with soil productivity and nutrient cycling, as well as many other ecosystem properties and services (24). Soil bacteria and fungi may share common resources, and competition for a substrate might induce the antagonism between bacteria and fungi (25). In addition, some soil-derived fungal and bacterial species may synthesize antibiotics (26,27), substantially affecting the species interactions. Bacteria may exhibit antifungal activity via producing volatile compounds, which were reported to inhibit the germination of fungal spores as well as hyphal growth and to change fungal morphology, enzyme activity, and gene expression (28). Species in highly competitive communities often grow less efficiently due to intense competition for shared resources (16,20). Antagonistic interactions of Pseudomonas fluorescens communities measured in vitro were related to bacterial root colonization and host plant protection, suggesting that increasing antagonistic interactions may cause negative BEF relationships (20). In turn, coexisting species resulting from niche partitioning by distinct resources may positively interact, which can increase functional community performance (14). For example, soil fungi may decompose the recalcitrant organic matter, e.g., cellulose and lignin, and bacteria may symbiotically utilize fungus-derived substrates (17). Since complex soil processes are driven by interactions among soil bacteria and fungi, etc., revealing the intrinsic linkages between microbial interactions and community functioning may facilitate the management of microbial communities for improving ecosystem service provisioning (29)(30)(31).
Here, we aimed to evaluate the potential factors which affect the strength of biodiversity-function relationships in complex terrestrial ecosystems. To address this issue, we conducted a large-scale soil survey ranging over different habitats, including an agricultural field, forest, wetland, grassland, and desert, along the Hexi Corridor (transect intervals of 1,257.6 km) (see Fig. S1 in the supplemental material), which is a representative of an oasis-desert ecotone in the arid regions of northwest China (32). This specific ecological environment of oases scattered along the narrow desert belt contains various ecosystems (32). We also took advantage of the strong changes in the soil biodiversity and processes that occur vertically along the soil profile. In total, 251 soil samples at 126 sites with two soil depth layers were selected in an agricultural field, forest, wetland, grassland, and desert. This cross-habitat environmental gradient offers a model system with an ecosystem and biodiversity gradient for the investigation of the relationships between biodiversity and ecosystem function (multinutrient cycling, for example). To control the effect of spatial scale, the sampling sites for each habitat were evenly distributed along the transect of the Hexi Corridor. We hypothesized that (i) microbial diversity-ecosystem function relationships would exhibit habitat-specific patterns and would be influenced by the community assembly processes and the balance between negative and positive species associations and that (ii) more negative than positive associations between bacteria and fungi would decrease the diversity-function link, along with decreasing neutral assembly processes. Our results may help predict and regulate biodiversity-driven ecosystem functioning and further develop proper land use strategies for improving the provision of key ecosystem services.

RESULTS
A total of 15,429,528 and 19,350,877 high-quality bacterial and fungal sequences were acquired from 251 cross-habitat soil samples along the Hexi Corridor (transect intervals of 1,257.6 km), which were, respectively, grouped into 25,981 and 21,698 operational taxonomic units (OTUs). Our results showed that the microbial a-diversity and multinutrient cycling index (MNC) in desert soils were lower than in other habitats (Fig. 1). The microbial a-diversities did not significantly differ between surface (depth, 0 to 15 cm) and subsurface (depth, 15 to 30 cm) soils in any of the habitats (see Fig. S2 in the supplemental material). In agricultural and forest soils, the MNC values were significantly higher in the surface soils than in the subsurface soils (Fig. S2). Moreover, there were significant differences in the microbial community compositions among the five habitats and between the two soil layers (Fig. S3A). In addition, a significant correlation between microbial a-diversities and community compositions was observed (Fig. S3B).
The relationships between microbial diversity and the MNC were explored in the five habitats and the two soil layers. There were significant positive relationships between microbial a-diversity and MNC only in agricultural and forest soils ( Fig. 2A). The microbial community composition was significantly correlated with the MNC in agricultural, forest, and wetland soils ( Fig. 2A). Similar trends were observed in the surface and subsurface layers in different habitats ( Fig. S4A and B). Concerning each component of multinutrient cycling, soil microbial diversity strongly correlated with most individual variables measured (Fig. 3). More negative correlations between microbial a-diversity and nutrient variables were observed in the grassland (n = 3) and desert (n = 4) than in other habitats (n # 2). The microbial community composition showed stronger correlations with most nutrient variables measured (Fig. 3) and higher correlation coefficients with MNC in 73% of habitat and soil layer combinations than a-diversities ( Fig. 2A and Fig. S4A), indicating the major role of microbial community composition in soil nutrient cycling. This observation was further confirmed by the multiple-regression model and variation partitioning analysis (Table S1).
We then estimated the microbial community assembly processes among the five habitats and between the two soil layers. The neutral community model explained a larger fraction of microbial community variation in agricultural soils (R 2 = 0.752) than in other habitats (R 2 , 0.70) (Fig. 2B). In addition, the degree of fit in the neutral community model was higher in surface layers than in subsurface layers of all habitats (Fig. S4C), suggesting that microbial communities were more governed by neutral processes in the surface layer, irrespective of habitat.
We further evaluated how the neutral assembly processes influenced the microbial diversity-MNC relationships. We observed significant relationships between the fit of the neutral model (R 2 ) and the correlation coefficients of the microbial diversity-MNC relationships (Fig. 4A), indicating that the habitat more governed by neutral processes exhibited a stronger diversity-MNC relationship. Second, we estimated how the balance between positive and negative bacterial-fungal associations correlated with diversity-MNC relationships. Interestingly, the proportions of negative associations between bacterial and fungal taxa were strongly and negatively (or positively) correlated with the correlation coefficients of the microbial diversity-MNC relationships ( Fig. 4B and Table 1). In other words, negative bacterial-fungal associations relative to positive ones modified the diversity-function link, whereas there were no significant relationships (P . 0.1) between the microbial diversity-MNC correlation coefficients and the negative associations among all taxa (Fig. S5A) within bacterial taxa (Fig. S5B) or fungal taxa (Fig. S5C). In addition, no significant correlations (P . 0.05) were detected between the soil pH-moisture ( Fig. S5D and E) and the microbial diversity-MNC relationships.
We applied a random-forest (RF) analysis to identify the major contributors to the microbial diversity-MNC relationships. We observed that the negative associations between bacterial and fungal taxa made a major contribution to predicting the microbial diversity-MNC relationships (Fig. 5). We then conducted structural equation modeling (SEM) to verify this observation (Fig. S6). The structural equation models had a FIG 1 Comparison of microbial a-diversities and soil multinutrient cycling among different habitats. The microbial a-diversity was calculated as the average value of bacterial and fungal diversities after minimummaximum normalization. The soil multinutrient cycling index (MNC) was calculated as the average value of soil organic carbon, dissolved organic carbon, microbial biomass carbon, nitrate-nitrogen, ammonium-nitrogen, microbial biomass nitrogen, available phosphorus, and available potassium contents after minimum-maximum normalization. The overall differences among habitats were estimated based on parametric one-way analysis of variance (ANOVA). In addition, different lowercase letters within panels indicate significant differences among the habitats (P , 0.05), which were found by multiple-comparison test after Kruskal-Wallis analysis. Error bars represent the standard errors. F, F-statistic.
FIG 2 Assessment of microbial diversity-multinutrient cycling relationships and community assembly processes in agricultural field, forest, wetland, grassland, and desert soil samples from the Hexi Corridor in northwest China. (A) Relationships between microbial diversity and multinutrient cycling. The microbial a-diversity was calculated as the average value of bacterial and fungal diversity after minimum-maximum normalization. The microbial community composition (b-diversity) was estimated using the first axis of the nonmetric multidimensional scaling analysis by combining the bacterial and fungal communities. Solid and dashed lines, respectively, denote the significant (P , 0.05) and nonsignificant (P . 0.05) Pearson's correlations. (B) Fit of Sloan's neutral model for analysis of microbial community assembly. The analysis was based on combining the bacterial and fungal communities. The solid blue line represents the best-fitting neutral model. The dashed line represents the 95% confidence intervals (CIs) around the best-fitting neutral model. OTUs within the CIs (black points) follow the neutral process. OTUs that occur more frequently than predicted by the model are shown in red, whereas those that occur less frequently than predicted are shown in blue. m indicates the estimated migration rate, and R 2 indicates the fit to the neutral model.
Drivers of the Diversity-Function Relationship good fit with the x 2 test, the root mean square error of approximation (RMSEA), the comparative fit index (CFI), and their P values. Overall, the models explained 66.4 and 67.8% of the variance found in the microbial diversity-MNC relationships for a-diversity and community composition, respectively ( Fig. 5C and D). In support of our earlier observations, we found the strongest and direct negative correlations between negative bacterial-fungal associations and the microbial diversity-MNC relationships. In addition, there was a negative relationship between the neutral processes and negative bacterial-fungal associations.

DISCUSSION
Soil microbial diversity promotes multifunctionality in natural terrestrial ecosystems (33, 34) but may be context dependent across complex ecosystems, with community assembly processes and balance between positive and negative interaction potentially controlling this context dependency. Therefore, characterizing determining factors is crucial for advancing the prediction and regulation of biodiversity-ecosystem function (BEF) relationships. Here, we reveal that community assembly processes and the balance between positive and negative bacterial-fungal associations were clearly linked with the microbial diversity-MNC relationships. The habitats less governed by neutral processes and dominated by negative bacterial-fungal associations exhibited stronger diversity-MNC relationships.
High soil microbial diversity can promote ecosystem multifunctionality and regulate multifunctionality resistance to climate change and fertilization in natural terrestrial ecosystems (33,35,36). Growing evidence for a strong link between soil biodiversity and multiple ecosystem functions has been reported (33,34,(36)(37)(38); however, these relationships have not yet been compared across habitats. Here, we observed that the microbial a-diversity and MNCs in desert soils were lower than in other habitats along the Hexi Corridor. This is supported by a previous study demonstrating that microbial communities in desert soils had the lowest levels of taxonomic diversity and gene abundances associated with nitrogen, potassium, and sulfur metabolism, compared to those of forests, grasslands, and tundra habitats (39). Moreover, the significant The shading from white to red represents low-to-high positive correlation, while the shading from white to blue represents low-to-high negative correlation. *, P , 0.05; **, P , 0.01; ***, P , 0.001. DOC, dissolved organic carbon; MBN, microbial biomass nitrogen; MBC, microbial biomass carbon; NH4, ammoniumnitrogen; NO3, nitrate-nitrogen; AP, available phosphorus; AK, available potassium; SOC, soil organic carbon. differences in the microbial community compositions among the five habitats and between the two soil layers were in agreement with recent studies that demonstrated strong habitat-specific patterns of microbial b-diversity in soil ecosystems (40,41). Specifically, we did not observe significant microbial diversity-MNC relationships in desert and grassland soils. Since oases were scattered along the narrow desert belt in arid regions, most of the grassland soil samples were collected from desert grassland in this study (40). In this case, the nonsignificant microbial diversity-MNC relationships in desert and grassland soils imply that desertification or grassland degradation might weaken the diversity-function relationship.
Positive effects of local-scale biodiversity on ecosystem functions have been demonstrated in theoretical, experimental, and observational studies across different types of ecosystems and habitats (21,24,34,36,38,42). Recently, the viewpoint that b-diversity is important in the context of multifunctionality was proposed (1). In the present study, we observed that community composition showed higher correlation coefficients with MNC than with a-diversity, suggesting the major role of microbial community composition on soil nutrient cycling. This supports the notion that community composition might be more important than community richness (43,44). In addition, the role of microbial community composition in regulating ecosystem functions is related to variations in local-scale diversity, which may scale up to large-scale changes in the provisioning of multiple ecosystem functions (45,46). Multifunctionality resistance to climate change and nitrogen fertilization is regulated by soil bacterial and fungal community compositions in natural ecosystems (47). Bacterial and archaeal b-diversities were strongly related to soil multinutrient cycling in a 30-year chronosequence of a reforestation ecosystem (48). S. Jiao et al. (41) demonstrated that soil multinutrient cycling in waterlogged rice fields was associated with bacterial b-diversity, Drivers of the Diversity-Function Relationship potentially attributable to the metabolic cooperation via syntrophy between bacterial groups under oxygen-limited conditions (41). Based on these cases, our results therefore support the perspective that community composition is closely linked to the provisioning of multiple ecosystem functions (multinutrient cycling, for example) (1). This may improve our understanding of BEF relationships under various natural and anthropogenic influences.
Biodiversity is a central topic in ecology, because the dramatic loss in biodiversity may reduce ecosystem functions and services (49). It is crucial to reveal the factors used to determine the strength of microbial diversity-function relationships in complex terrestrial ecosystems, considering the intrinsic linkages between assembly processes and species interactions. Understanding the assembly mechanisms of belowground microbial communities is crucial for understanding the maintenance and generation of terrestrial microbial diversity (50)(51)(52). Previous studies showed that microbial biogeography exhibited strong habitat-specific patterns (40,41). In the present study, we found that the assembly of microbial communities in agricultural soils was more governed by neutral (e.g., stochastic) processes than that in other natural habitats, indicating that long-term cultivation and human-managed activities (53) might enhance the stochastic influx and dispersal of microorganisms. Our result supported prior studies reporting that stochastic processes were stronger in crop fields than in grassland (54). Additionally, we observed that neutral processes were stronger in surface soils than in subsurface layers, irrespective of habitat. This might be due to the difference in dispersal between surface and subsurface soils (55). Particularly, available nutrient substrates usually decrease with soil depth because of the reduced input of plant litter and root exudates (56). High available nutrients introduce more resources to soil microbes, which can improve the ability of microorganisms to disperse and in turn increase the dominance of stochastic (e.g., neutral) processes (3). Therefore, our results suggest that soil depth may affect the microbial community assembly processes.
Species interactions are considered to be deterministic (niche-based) assembly processes governing community structure (6,57,58). In the present study, the SEM results showed a negative relationship between the neutral processes and negative bacterial-fungal associations, indicating that negative relationships tended to be lower when communities were driven primarily by neutral processes. That is, positive bacterial-fungal associations were more widespread under stronger stochastic processes. This may be supported by a previous study demonstrating that soil microbial cooccurrence associations are higher when communities are driven primarily by stochastic processes (e.g., dispersal limitation) (59). Stochastic processes involve random birth, death, and dispersal events; therefore, species tend to cooccur when random changes and increased influxes of species are not associated with environmentally derived fitness (6, 7). Previous studies have demonstrated that community assembly processes were associated with the microbial interactions (e.g., cooccurrence associations based on the correlation network analysis) in soil (59) and river (60) systems. In addition, a global soil microbiome study provided evidence for strong bacterial-fungal antagonism, suggesting the role of species interactions in shaping microbial communities (12). On this basis, our findings more specifically uncover the intrinsic linkages between assembly processes and microbial associations (e.g., positive or negative) and suggest that the balance between positive and negative bacterial-fungal associations is strongly related to community assembly processes.

Drivers of the Diversity-Function Relationship
Furthermore, we observed that the neutral processes substantially influenced the magnitude and direction of the diversity-MNC relationships. The habitats more governed by neutral processes exhibited stronger positive or negative diversity-MNC relationships for a-diversity or community composition, respectively. This may be explained by two mechanisms: (i) immigration can add species or change species composition in a way that increases the biodiversity effect on functions through the sampling effect if immigration from the regional pool brings new species carrying traits that affect ecosystem functioning and that were not present in the initial community (10) and (ii) a higher relative balance of a deterministic process (and thus the lower importance of neutral ones) can decrease the biodiversity effect on ecosystem function by a dilution effect if species selected by deterministic processes are not the ones affecting ecosystem functioning (10). Thus, a lower influence of deterministic processes (and thus a higher influence of stochastic ones) might reduce this dilution effect and strengthen the BEF relationship, as observed in the present study. A previous study used a model to assess how the relative balance between stochastic and deterministic processes affect a generic biogeochemical function and revealed that higher dispersal led to decreases in biogeochemical function due to the increased abundance of poorly adapted organisms (61). The inconsistency might be attributed to the distinct models, and we considered only the effect of microbial diversity on biogeochemical functions. Uniquely, our study builds the linkage between the mechanisms underlying community assembly and the BEF relationships and suggests that neutral processes governing the distinct assembly patterns of a cross-habitat microbial community may enhance the functional contributions of ecological communities. This highlights the potential roles of community assembly mechanisms in generating and sustaining the multiplenutrient cycling of terrestrial ecosystems.
Species interactions play important roles in stimulating ecosystem processes (62). In the present study, we observed that the proportions of negative associations between bacterial and fungal taxa were clearly linked with the microbial diversity-MNC relationship when we considered environmental factors and community assembly processes. It indicated that the balance between positive and negative bacterial-fungal associations underpinned the context dependency of the microbial diversity-MNC relationship. Negative species associations may be due to antagonistic biological interactions (63), including competition. It is thus conceivable that a high number of antagonistic or competitive interactions between bacterial and fungal taxa led to the weak microbial diversity-MNC relationships observed in grassland and desert habitats, which was supported by the SEM results. Competition was a major type of interaction between fungi and bacteria in soil (19,20,64,65). For example, some soil-derived fungal and bacterial species may synthesize antibiotics and showed antagonistic effects (26,27) which had consequences for microbial community assembly (12). A study has shown that strongly hierarchical competitive network communities comprising strong competitors exhibit a negative diversity-function relationship (19). Positively interacting species can increase functional community performance due to niche partitioning by distinct resources (15,16). Soil organisms with similar environmental preferences may form strongly connected ecological clusters in ecological networks, with major implications for ecosystem functioning (66,67). Previous studies showed that fungal-bacterial interkingdom associations may enhance ecosystem functioning related to nutrient cycling in grasslands (31) and promote plant health in the model plant Arabidopsis (68). Thus, these studies may support our conclusion that the balance between positive and negative associations were connected to the link between soil microbial diversity and multinutrient cycling.
One potential limitation of this study is that the negative/positive associations between taxa should not be interpreted as a proof of competition/facilitation, which was based only on correlation (69). Correlation analyses are only a simplistic representation of a complex system, although they are frequently used to investigate microbial interconnection patterns (70)(71)(72). In addition, species associations that are based on correlations can yield spurious results and cannot be automatically interpreted as interactions (63). Consequently, it may not be possible to comprehensively depict the microbial interactions under real-world conditions. However, the information about negative/positive correlations between taxa is still essential for estimating potential species interrelationships within complex environments and, in turn, for revealing the influence of microbial interconnection complexity on biodiversity-driven ecosystem functioning (67).
Conclusions. Our study highlights that the balance between positive and negative bacterial-fungal associations is clearly linked with the strength of the relationships between soil microbial diversity and multiple nutrients cycling across different habitats. Changes in both the level of neutrality and the proportion of negative bacterialfungal associations are linked to the magnitude and the direction of the diversity-MNC relationship. These findings reveal the potential factors underpinning the context dependency of the microbial diversity-MNC relationship, which makes it possible to predict the ecological consequences for the biodiversity-function relationship in the future by uncovering how soil microbial assembly and associations are likely to respond to the climate and land use changes. For example, the nonsignificant microbial diversity-MNC relationship in desert soils suggests that desertification might weaken the microbial diversity-function relationship, along with increasing negative bacterial-fungal associations. The strongest neutral process of microbial assembly in agricultural soils implies that the conversion from natural habitats to agricultural fields (e.g., the conversion of grassland for agricultural crops, the conversion of forests for crop and wetland degradation) might enhance the stochastic influx and dispersal of microorganisms and strengthen the biodiversity-MNC relationship. Overall, our results represent a considerable advancement in facilitating the management of the functioning of microbial communities for improving human well-being, in addition to informing strategies based on community assembly and microbial interactions.

MATERIALS AND METHODS
Site description. The study site extended from 36°569N to 40°349N and from 94°379E to 103°319E along the Hexi Corridor in the northwestern portion of Gansu Province and to the west of the Yellow River in China (see Fig. S1 in the supplemental material). The Hexi Corridor is a long belt between the South Mountains (including Mt. Qilian and Mt. Aerjin) and the North Mountains (including Mt. Mazong, Mt. Heli, and Mt. Longshou). Due to the distribution of oases scattered along the narrow desert belt, the Hexi corridor contains a variety of soil ecosystems (32).
Here, we selected five habitats: an agricultural field, forest, wetland, grassland, and desert, according to a vegetation map at a scale of 1:1,000,000 (Data Centre for Resources and Environmental Sciences, Chinese Academy of Sciences [RESDC], http://www.resdc.cn). The dominant species in these habitats included Zea mays (agricultural field), Calligonum spp., Stipa spp., Leymus spp., and Achnatherum spp. (wetland, grassland, and desert), and Populus spp. (forest). The dominant soil types were aripsamment and calciorthids, which have a loose structure and low organic matter content. The climate of this region is predominantly semiarid.
Sample collection. Field sampling was conducted during July and August 2017 near the period of the highest aboveground plant biomass. To ensure appropriate spatial scale, the sampling sites for each habitat were evenly distributed along the transect of the corresponding region. This meant that some sampling sites for different habitats were spaced less than 5 km apart; hence, some of these sites overlapped on the map (Fig. S1). In total, 126 sites were selected; 37 were from an agricultural field, 28 were from a forest, 15 were from a wetland, 26 were from a grassland, and 20 were from a desert.
At each site, three plots were sampled; each plot had an area of 100 m 2 . Five soil cores (2.5-cm diameter) were combined per plot and were taken at depths of 0 to 15 cm and 15 to 30 cm (56,73). The soil cores from the three plots for a given soil depth layer were then mixed thoroughly to generate the final soil samples. All soil samples were delivered to the laboratory in sterile plastic bags on dry ice and were sieved through a 2.0-mm mesh to remove plant debris and rocks. A portion of each soil sample was stored at 4°C for the analysis of nutrient factors. Aliquots of soil samples were stored at 220°C for subsequent DNA extraction. One subsurface desert sample was abandoned due to DNA extraction failure. Therefore, a total of 251 soil samples at 126 sites with two soil depth layers were used for this study. Standard testing methods were applied to measure soil pH, moisture, soil organic carbon, dissolved organic carbon, microbial biomass carbon, nitrate-nitrogen, ammonium-nitrogen, microbial biomass nitrogen, available phosphorus, and available potassium, as previously described (50,71,74). Detailed descriptions are provided in the supplemental material. We acknowledge that many important environmental parameters (e.g., the types of carbon/nitrogen/phosphorus, the hydrology, temperature, oxygen stress, microsite structure, anion/cation/micronutrients, etc.) were not measured in this study, and these should be considered in future work.

Drivers of the Diversity-Function Relationship
Microbial DNA processing. Bacterial and fungal diversity scores were obtained via high-throughput sequencing of the PCR amplicons of the 16S rRNA and internal transcribed spacer (ITS) genes. Briefly, the total genomic DNA was extracted from the soil samples (0.5 g) using a FastDNA SPIN kit for soil (MP Biochemicals, Solon, OH, USA). The microbial communities were profiled by targeting the V4-V5 region of the 16S rRNA gene for bacteria and the ITS1 region of the 18S rRNA gene for fungi. The target sequences were amplified by PCR using the primer pairs 515F/907R (bacteria) and ITS5-1737F/ITS2-2043R (fungi) (75,76). Sequencing was conducted on an Illumina HiSeq2500 platform (Illumina Inc., San Diego, CA, USA).
The acquired sequences were filtered for quality according to the method of J. G. Caporaso et al. (77). Sequences were assigned to their corresponding samples according to the barcode and then quality trimmed with a threshold of average Phred quality scores of higher than 20. Any chimeric sequences were removed with the USEARCH tool based on the UCHIME algorithm (78). The sequences were split into groups according to their identity and assigned to operational taxonomic units (OTUs) at a 3% dissimilarity level using the UPARSE pipeline (78). The OTUs with no more than two sequences were removed, and their representative sequences were classified within the SILVA database (release 128) for bacteria (79) and UNITE plus INSD (UNITE and the International Nucleotide Sequence Databases; release 7) for fungi (80). Counts of individual OTUs were scaled by the total number of reads in each sample to account for sequencing biases using the R package DESeq2 (81). This measure of normalized abundance allows samples with various read counts to be compared (82) and is widely applied to high-throughput sequencing data (82,83).
Soil microbial diversity analysis. To obtain a multidiversity index, we combined soil microbial diversity characteristics by averaging the standardized scores of bacterial Shannon diversity and fungal richness. The scores standardized to a common scale ranging from 0 to 1 were calculated according to the following formula: STD = (X 2 X min )/(X max 2 X min ), where STD is the standardized variable and X, X min , and X max are the target variable and its minimum and maximum values across all samples, respectively. This multidiversity index is largely used and accepted in the current biodiversity function literature (34,37,84). The microbial b-diversity was quantified using the first axis of a nonmetric multidimensional scaling (NMDS) analysis of Bray-Curtis dissimilarities by combining the bacterial and fungal communities (85). Here, the analysis was based on combining the bacterial and fungal communities, including all of the bacterial and fungal OTUs. Before the combination, the bacterial and fungal communities were standardized to a total abundance of 1.
Evaluation of soil ecosystem function. Ecosystems perform multiple simultaneous functions and services, rather than a single measurable process. Given that nutrient cycling is the most important soil ecosystem process for supporting human welfare (62,86), we estimated the multinutrient cycling index (MNC) to evaluate soil ecosystem functions (41,48). This index comprised information for 8 soil nutrient variables in relation to carbon (soil organic carbon, dissolved organic carbon, and microbial biomass carbon), nitrogen (nitrate-nitrogen, ammonium-nitrogen, and microbial biomass nitrogen), phosphorus (available phosphorus), and potassium (available potassium) cycling. These variables constitute an integrated proxy for nutrient cycling and are important determinants of ecosystem functioning in terrestrial ecosystems (33,34,84). For example, nitrogen and phosphorus are the nutrients that most frequently limit primary production in terrestrial ecosystems (87). In addition, potassium is the third essential macronutrient required by plants; it participates in a range of biological activities, such as protein synthesis, enzyme activation, and photosynthesis, that maintain or improve plant growth (88). We acknowledge that some important functions, such as soil process rates, are inevitably unmeasured, and future studies are encouraged to include more essential functions for comprehensive understanding of ecosystem functioning. To derive a quantitative MNC value for each site, we averaged the standardized scores (a common scale ranging from 0 to 1) of all individual nutrient variables. This method was used to quantify soil multinutrient cycling because it is a straightforward and interpretable measure of a community's ability to sustain multiple functions simultaneously (33,34,84).
Ecological analysis. The NMDS analysis was performed to visualize the sample relationships across different habitats. An analysis of similarities (ANOSIM) was used to determine significant differences in microbial community composition across different habitats, performed using the anosim function in the vegan package in R (89). Multiple-comparison post hoc tests for Kruskal-Wallis analysis were used to test for significant differences in microbial diversity and soil multinutrient cycling among different habitats, performed using the kruskal function in the agricolae package in R (90). Pearson's correlation analysis was used to estimate the relationship between microbial diversity and multinutrient cycling, performed using the cor.test function in the stats package in R (91).
Neutral modeling. A Sloan neutral community model was used to determine the contribution of neutral processes to microbial community assembly (92). The model predicts that abundant taxa are more likely to be dispersed by chance and widespread across a metacommunity, while rare taxa are lost in different local communities due to ecological drift. The neutral model is fit to the relationship between the frequency with which taxa occur in a set of local communities and their abundance across the wider metacommunity by estimating the parameter describing the migration rate (m), a measure of dispersal limitation. Higher m values indicate that microbial communities are less dispersal limited (92,93). The formula (59) is Freq ¼ 1 2 I 1=N; N Â m Â p; N Â mð1 2 pÞ ½ , where Freq is taxon occurrence frequency, N is the number of individuals per community, p is the taxon relative abundance, and I[] is the probability density function of the beta distribution. R 2 indicates the fit of the parameter based on nonlinear least-squares fitting. The overall fit of the model to the observed data was assessed by comparing the sums of squares of residuals, SS err , with the total sum of squares, SS total : model fit = 1 2 SS err /SS total (generalized R-squared) (93). Higher R 2 values indicate a higher contribution of neutral processes to microbial community assembly. In the present study, we used the fit of the neutral model (R 2 ) to infer the neutral assembly processes. One point should be noticed, i.e., that stochastic processes do not exactly incorporate a neutral process, although a few recent researchers have applied neutral-theory-based process models to infer the stochastic processes (57,94). Here, the analysis was based on combining the bacterial and fungal communities, which were, respectively, standardized to a total abundance of 1.
Microbial association analysis. To explore the potential interactions among species, we estimated the associations among all the members of the bacterial and fungal communities in different habitats and soil layers. Robust correlations were estimated via Spearman's correlation analysis with false-discovery rate (FDR)-corrected P values of ,0.01, which were used to reflect the negative (Spearman's correlation coefficient [r ] , 0) or positive (r . 0) associations among microbial taxa. The proportion of negative associations meant negative associations divided by the total associations. To avoid random effects of rare taxa, only taxa detected in more than 60% of the soil samples in each habitat (e.g., different habitats and layers) were used for the correlation analysis (95). To test whether the outcomes were sensitive to the choice, we also estimated the microbial associations based on the 50% threshold. Similar results with the 60% threshold were observed (data not shown), indicating that the outcomes were not sensitive to the choice of threshold.
RF modeling. We first evaluated Pearson's correlations between the strength of microbial diversity-MNC relationships (correlation coefficients) and (i) the proportions of negative associations between bacterial and fungal taxa, within all taxa, within bacterial taxa, and within fungal taxa; (ii) the neutral community assembly processes (R 2 of the neutral model); and (iii) environmental variables, including soil pH and moisture. Additionally, random-forest (RF) analysis was performed to identify the main factors influencing the microbial diversity-MNC relationships (33,96). In the RF models, negative associations between bacterial and fungal taxa, neutral community assembly processes, soil pHs, and moisture levels served as predictors for the correlation coefficients of the microbial diversity-MNC relationships. To estimate the importance of these variables, we used percentage increases in the mean squared error (MSE) of variables: higher MSE percentages imply more important variables (97). The significance of the model was assessed with 5,000 permutations of the response variable by using the A3 package (98). Similarly, the significance of each predictor on the response variables was assessed with 5,000 trees by using the rfPermute package (99).
SEM. We then used structural equation modeling (SEM) to evaluate the direct and indirect effects of different factors on the strengths of microbial diversity-MNC relationships. The first step in SEM requires establishing an a priori model based on the known effects and relationships among the drivers (Fig. S6). The negative associations between bacterial and fungal taxa, neutral community assembly processes, soil pH, and moisture were considered in the model. Since our sampling sites for each habitat were evenly distributed along the transect of the Hexi Corridor to control the effect of spatial scale, the spatial and climatic variables were not included in the model. We fitted the full model containing all potential paths of our a priori model (Fig. S6) and then simplified the model by removing the variable (e.g., soil pH) without any significant relationship (44). Each path removal was accepted if the model quality-based Akaike information criterion (AIC) was improved. The goodness of fit of structural equation models was checked using the following: the x 2 test, the root mean square error of approximation (RMSEA), and the comparative fit index (CFI). The model has a good fit when the CFI value is close to 1, RMSEA values are closer to 0, and x 2 and RMSEA P values are high (traditionally .0.05) (100). With a good model fit, we were free to interpret the path coefficients of the model and their associated P values. A path coefficient is analogous to the partial correlation coefficient and describes the strength and sign of the relationship between two variables. SEM was conducted with the lavaan package (101).
Data accessibility. The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (107) and in the Beijing Institute of Genomics (BIG) Data Center (108), BIG, Chinese Academy of Sciences, under BioProject accession no. PRJCA004036 and are publicly accessible at http://bigd.big.ac.cn/gsa.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, DOCX file, 0.04 MB.