Distinctive pattern and mechanism of precipitation changes affecting soil microbial assemblages in the Eurasian steppe.

Precipitation may increase or decrease by different intensities, but the pattern and mechanism of soil microbial community assembly under various precipitation changes remain relatively underexplored. Here, although ±30% precipitation caused a small decrease (∼19%) in the within-treatment taxonomic compositional dissimilarity through the deterministic competitive exclusion process in a steppe ecosystem, ±60% precipitation caused a large increase (∼35%) in the dissimilarity through the stochastic ecological drift process (random birth/death), which was in contrast with the traditional thought that increasing the magnitude of environmental changes (e.g., from +30% to +60%) would elevate the importance of deterministic relative to stochastic processes. The increased taxonomic dissimilarity/stochasticity under ±60% precipitation translated into functional dissimilarity/stochasticity at the gene, protein, and enzyme levels. Overall, our results revealed the distinctive pattern and mechanism of precipitation changes affecting soil microbial community assembly and demonstrated the need to integrate microbial taxonomic information to better predict their functional responses to precipitation changes.

Introduction

Soil microbial communities harbor perhaps the highest biodiversity on the planet, and they are responsible for many important ecosystem functions such as carbon (C) and nitrogen (N) cycling (Fierer et al., 2012; Torsvik et al., 2002). Anthropogenic climate change may alter microbial community diversity and composition and their ecosystem functions, which can result in positive or negative feedbacks to climate change (Guo et al., 2020; Heimann and Reichstein, 2008; Singh et al., 2010). Although changes in both temperature and precipitation are key components of climate change, the effects of precipitation changes on soil microbial communities remain underexplored, especially compared with the large number of studies that have focused on the effects of temperature increases (Karhu et al., 2014; Nie et al., 2012; Ning et al., 2020). Precipitation may increase or decrease by different intensities as a consequence of climate change, but few studies have investigated the patterns and mechanisms underlying the response of soil microbial communities to these precipitation changes. Because water makes up more than 60% of the weight of cellular organisms, precipitation should be the most important ecological factor structuring terrestrial biological communities (Jorgensen, 2009). Consistently, precipitation has been shown to be more important than other ecological factors (e.g., soil nutrient content) in structuring aboveground plant communities over large spatiotemporal scales (Bai et al., 2004, 2008). Precipitation has also been found to play some roles in structuring belowground microbial communities across some natural and experimental precipitation gradients (Averill et al., 2016; Delgado-Baquerizo et al., 2018; Hawkes et al., 2011); however, soil pH rather than precipitation (or soil water content) is often found to be the most important factor in structuring soil bacterial communities over different spatiotemporal scales (Fierer, 2017), and various types of anthropogenic environmental changes are found to affect soil microbial communities primarily through mediating soil pH rather than soil water content (Zhang et al., 2014; Zhou et al., 2020a). Therefore, the direct effect and critical role of precipitation on soil microbial diversity might be underestimated in previous studies. The underestimation of the importance of precipitation in structuring soil microbial communities suggests that precipitation might affect soil microbial assemblages via distinctive pathways and mechanisms that have yet to be discovered.

Ecologists have long been interested in the mechanisms of species coexistence within a community (Chesson 2000; Hutchinson, 1959). The traditional niche theories propose that different species coexist in a community because they have different niches, e.g. they are limited by different resources (Harpole and Tilman, 2007; Pacala and Tilman, 1993). Within a community, different species have different fitness, and their abundances are consistent with their fitness. As environmental condition changes, their fitness will change inconsistently, leading to the alteration of their abundances and the community structure. Generally, the niche theories focus on the deterministic effects of selection. In recent years, ecologists have proposed ecological neutral theories, which posit that changes in community structure are not related to fitness differences among species but are instead affected more by stochastic processes associated with ecological drift, migration, and speciation (Hubbell 1979, 2001). Both deterministic processes associated with differences in fitness among species, such as competitive exclusion and environmental filtering, and stochastic processes that are independent of differences in fitness among species, such as ecological drift and random dispersal/colonization, have been found to be potential drivers of microbial community assembly (Hubbell, 2001; Tilman, 2004). Anthropogenic environmental changes can affect community assembly by mediating both deterministic and stochastic processes (Zhang et al., 2016; Zhou and Ning, 2017). For example, soil acidification has been shown to drive soil microbial community assembly primarily through deterministic processes, whereas nutrient enrichment drives the assembly of both producer and animal communities in freshwater ponds primarily through stochastic processes (Chase, 2010; Zhang et al., 2011). Increases in the intensity of environmental change are thought to impose strong selection pressures on communities and favor taxa with higher fitness, thus increasing the relative contribution of deterministic processes compared with stochastic processes (Zhang et al., 2011). Given that precipitation might be the most important factor affecting soil microbial assemblages, we propose a hypothesis about the distinctive pattern and mechanism by which precipitation changes driving microbial community assembly (Figure 1). When the amount of precipitation increases (or decreases) below a certain threshold and the carrying capacity (the potential maximum community size) of the soil habitat changes by only a small extent (Figure 1A), the importance of deterministic processes (namely, competitive exclusion) in soil microbial community assembly would increase (Figure 1B), and this would lead to increases in the abundance of, for example, hygrophilous (or drought-tolerant) microbial species. This would in turn result in a decrease in the within-treatment compositional dissimilarity (Figure 1C) and an increase in the predictability of community functions (Figure 1D). When the amount of precipitation increases (or decreases) above a certain threshold, the soil habitat can support more (or fewer) microbial individuals (Figure 1A), and this would stimulate the random birth (or death) process (Figure 1B). As a consequence, the within-treatment compositional dissimilarity and the unpredictability of community functions would increase (Figures 1C and 1D). In other words, the above-threshold changes in precipitation affect microbial community assembly by stimulating the stochastic process of ecological drift (random birth/death), and these stochastic changes in community composition might contribute to explaining why the importance of precipitation has been underestimated in previous studies.

Figure 1

Conceptual framework of precipitation changes affecting soil microbial assemblages

(A–D) Conceptual framework of how precipitation changes affect community size (A), the relative contributions of deterministic relative to stochastic processes in community assembly (B), and the within-treatment community compositional dissimilarity (C), leading to changes in the predictability of community function (D).

To test our hypothesis, we performed a 5-year precipitation experiment with five treatments (the ambient control, ±30% precipitation, and ±60% precipitation) in a steppe ecosystem in Inner Mongolia, China, that is floristically and ecologically representative of much of the Eurasian steppe region (Li, et al., 1988). According to the historical growing-season precipitation in the area, the amount of precipitation in the driest (or wettest) year was approximately 60% below (or above) the long-term average precipitation. Thus, ±60% precipitation was used for the above-threshold precipitation changes, and ±30% precipitation was used for the below-threshold changes. We quantified the taxonomic diversity of all soil prokaryotic communities (bacteria and archaea) using 16S rDNA deep amplicon sequencing and measured soil microbial functional diversity with metagenome and metaproteome sequencing. We used a null model method to compare the observed communities with the stochastically assembled communities (Chase, 2010; Zhou et al., 2014) and determine the relative contributions of deterministic and stochastic processes in driving community assembly. An increasing number of studies have shown that stochastic processes are key drivers of microbial community assembly, and many have attributed this pattern to the large number of rare species in microbial communities (Sloan et al., 2006). To assess the contribution of rare species to determining the importance of deterministic versus stochastic community assembly, we varied the number of rare operational taxonomic units (OTUs) included in the analysis of the deep amplicon data. To determine whether the stochastic process through which the above-threshold precipitation change affects microbial community assembly is a distinctive mechanism, we re-analyzed the soil microbial community composition data from a previous multifactorial experiment conducted in a similar steppe ecosystem in Northern China. In this experiment, a total of 13 environmental changes, including watering (above-threshold precipitation increase), warming, N addition, phosphorus addition, mowing, and some of their combinations, were examined. We determined whether the relative contributions of stochastic processes relative to deterministic processes in microbial community assembly under watering (and coupled with other environmental changes) were larger than those under the other types of environmental changes.

Results and discussion

±30% and ±60% precipitation treatments were used as the below- and above-threshold changes, respectively

There were five precipitation treatments (−60%, −30%, 0 (control, natural precipitation), +30%, and +60%), with four replicates for each treatment. Soil samples were collected after the 5-year experiment. Increases in precipitation significantly increased the soil water content, soil available N content, plant biomass, and microbial biomass C and N content and decreased soil temperature (p < 0.05; Figure 2); however, increases in precipitation did not significantly affect soil pH, soil organic C and total N content, dissolved organic C content, and plant richness. Although −30% precipitation reduced the soil water content by only 10.92%, −60% precipitation reduced the soil water content by 26.57% (>2 × 10.92%; Figure 2A); +30% precipitation increased the soil water content by only 4.44%, but +60% precipitation increased it by 15.46% (>2 × 4.44%; Figure 2A). Similarly, although −30% precipitation reduced the microbial biomass N content by only 8.82%, −60% precipitation reduced it by 55.45% (>2 × 8.82%; Figure 2K). These results indicated that the soil water content and microbial biomass appear to have a certain buffering capacity to the ±30% precipitation treatments. Thus, ±30% precipitation was used as the below-threshold precipitation change, which altered the carrying capacity (the potential maximum community size) of the soil habitat by only a small extent, and ±60% precipitation was used as the above-threshold precipitation change, which greatly altered the carrying capacity, supporting our hypothesis in Figure 1A.

Figure 2

Effect of precipitation changes on the biotic and abiotic indices

Below- and above-threshold precipitation changes decreased and increased the within-treatment compositional dissimilarity, respectively

To detect the relatively rare microbial species and accurately measure microbial community composition, deep amplicon sequencing targeting the 16S rDNA sequences was conducted for each of the four replicate communities of all five precipitation treatments. A total of 716,472–851,990 sequences were obtained from each of the 20 bacterial communities, and 148,106–244,669 sequences were obtained from each of the 20 archaeal communities. These sequences were classified into OTUs using a >97% similarity threshold. A total of 13,205 bacterial OTUs and 214 archaeal OTUs were observed across all 20 communities, with 9,420 ± 91 (mean ± standard error) bacterial OTUs and 96 ± 4 archaeal OTUs in each community. One-way ANOVA revealed that the experimental precipitation gradient had a nonsignificant effect on bacterial and archaeal OTU richness (p > 0.05), and this result was robust to increases in the number of rare OTUs included in the analysis (e.g., including OTUs with >1% or >0.001% relative abundance; Table S1). These results indicated that experimental precipitation did not affect soil microbial OTU richness.

Both weighted and unweighted Unifrac distances, which represented the phylogenetic distances between communities, were calculated for soil bacterial and archaeal communities (Lozupone et al., 2007). PERMANOVA (permutational multivariate analysis of variance) is a permutation method for testing the response of variables to factors in an ANOVA experimental design on the basis of any distance measure, and PERMDISP (permutational analysis of multivariate dispersions) is a method for comparing the multivariate dispersion among groups on the basis of any distance or dissimilarity measure. Both PERMANOVA and PERMDISP revealed that the precipitation changes did not significantly affect the phylogenetic structure of soil bacterial and archaeal communities (p > 0.05; Table S2). Bray-Curtis dissimilarity based on the relative abundances of OTUs was calculated to characterize variation in community composition. PERMANOVA revealed that the precipitation gradient did not result in any directional/systematic changes in bacterial or archaeal community composition (p > 0.05; Figure S1). However, PERMDISP revealed that the precipitation gradient had a significant effect on the within-treatment compositional variation (p < 0.01; Figure S1). Overall, these results indicated that the precipitation gradient had nondirectional/nonsystematic effects on microbial community composition. In other words, although the overall composition of all replicate communities of a treatment did not differ from that of the control, the composition of a single replicate community of a treatment might differ from the overall composition of the control, and the four replicate communities of a certain treatment diverged from the overall composition of the control in an inconsistent and unpredictable manner. This is an interesting observation because experimental treatments have often been observed to alter community composition directionally/systematically rather than nondirectionally/nonsystematically (Chase, 2007; Zhang et al., 2011, 2019). For example, the soil microbial community composition of N-enrichment plots was found to shift directionally rather than nondirectionally relative to control plots according to the ordination plots (Zhang et al., 2019).

To test whether the nonsystematic response of community composition to the precipitation gradient stemmed from the large number of rare OTUs, we analyzed the community data by varying the number of rare OTUs included in the analysis (e.g., including OTUs with >1%, >0.001%, and even lower relative abundances). The nonsignificant PERMANOVA results and the significant PERMDISP results were highly robust (Table S2), which suggested that the common and rare OTUs did not differ in their response to precipitation changes; this finding also confirmed that the response of community composition is not systematic.

Among all 13,205 bacterial OTUs, the relative abundances of 47.38% of the OTUs were significantly heterogeneous across the five precipitation treatments (Levene’s test: p < 0.05; Table S3), which is consistent with the significant PERMDISP results; the relative abundances of only 6.11% of the OTUs were significantly correlated with the soil water content (Pearson correlation: p < 0.05; Table S3), which is consistent with the nonsignificant PERMANOVA results. Similarly, among the 214 archaeal OTUs, 75.70% of the relative abundances of OTUs were significantly heterogeneous across the five precipitation treatments, and only 3.74% of the relative abundances of the OTUs were significantly correlated with the soil water content (Table S4). Overall, these results confirmed the nonsystematic response of microbial community composition to the precipitation gradient.

PERMDISP revealed a consistent pattern for both soil bacterial and archaeal communities in which ±30% precipitation resulted in a relatively small decrease in the within-treatment compositional dissimilarity compared with the control, whereas ±60% precipitation resulted in a relatively large increase in the within-treatment compositional dissimilarity (Figures 3A and 3B). This pattern was robust to increases in the number of rare OTUs included in the analyses (Figure S2). Consistent with the decrease in the within-treatment compositional dissimilarity under the ±30% precipitation treatments, 29.53% of bacterial OTUs and 63.55% of archaeal OTUs exhibited their lowest variance in relative abundance under the two treatments. Consistent with the increase in the within-treatment compositional dissimilarity under the ±60% precipitation treatments, 28.50% of bacterial OTUs and 41.59% of archaeal OTUs showed their largest variance in relative abundance under these two treatments (Tables S3 and S4). This pattern supported our hypothesis in Figure 1C that below-threshold precipitation changes (±30% precipitation) decreased the within-treatment compositional dissimilarity, whereas above-threshold precipitation changes (±60% precipitation) increased the within-treatment dissimilarity.

Figure 3

Precipitation changes affect community composition and ecological drivers

(A–D) Effect of precipitation changes on within-treatment compositional dissimilarity (A, B) and the SES (standard effect size (C, D) value of soil bacterial and archaeal communities revealed by PERMDISP.

For multiple comparisons, significant (p < 0.05) differences are differentiated by upper case letters. Error bars represent standard errors. OTUs with >0.01% relative abundance were included in the analysis. See also Figures S1–S5 and S8, Data S1, and Tables S1, S2, S3, S4, S5, S6, S7, and S8.

Below- and above-threshold precipitation changes increased the relative importance of deterministic and stochastic processes, respectively

A Mantel test was used to identify the factors (plant community indices and soil physicochemical indices) explaining variation in microbial community composition. Variation in bacteria and archaea community composition was significantly correlated with certain indices, including the soil water content, temperature, pH, available N content, and plant community composition (p < 0.05; Table S5). However, the r-value of the Mantel test was generally <0.5, meaning that most of the variation was not explained by these soil and plant indices. Consistently, a variation partitioning analysis indicated that substantial portions (75.7–85.7%) of the variation in bacterial and archaeal community composition could not be explained by the measured environmental variables (Figure S3), suggesting that stochastic processes and/or unmeasured environmental variables played an important role in the assembly of soil microbial communities (Ramette and Tiedje, 2007).

To analyze the relative contributions of deterministic versus stochastic processes in driving soil microbial community assembly in each treatment, we tested for significant differences between the observed within-treatment dissimilarity (β diversity) of the bacterial and archaeal communities and the expected dissimilarity of stochastically assembled communities using a null model method (Chase, 2010; Zhou et al., 2014). This method controlled for the effects of α and γ diversity, and thus permitted the effects of ecological processes on β diversity to be determined. No significant differences between the observed and expected dissimilarity were observed for control communities (p > 0.05; Table S6), suggesting that stochastic processes were the primary driver of the assembly of control microbial communities in this soil ecosystem. Microbial community assemblages have often been found to be governed by stochastic processes in other ecosystems (Sloan et al., 2006; Zhou and Ning, 2017), which has often been thought to stem from their large number of rare species. However, the nonsignificant differences between the observed and expected dissimilarity were robust to increases in the number of rare OTUs in the analyses (Table S6), indicating that the stochastic nature of microbial community assembly was not primarily caused by the large number of rare species. Thus, the stochastic nature of microbial community assembly is likely determined by other characteristics of microorganisms, such as their small body sizes and short generation times; however, this hypothesis requires confirmation.

For the other four precipitation treatments, the differences between the observed within-treatment dissimilarities and the expected dissimilarities of stochastically assembled communities were not significant (p > 0.05; Table S6) in 35 out of the 40 cases (2 microbial groups × 5 OTU rarity levels × 4 precipitation gradients), demonstrating that stochastic processes were the primary driver of soil microbial community assembly at different precipitation levels. Thus, differences were significant in only 5 out of 40 cases (p < 0.05; Table S6), suggesting that experimental precipitation changes only weakly affected the relative contributions of stochastic versus deterministic processes to driving soil microbial community assembly.

To further investigate the effects of precipitation changes, the community data of all plots of all treatments were analyzed together using the null model method to calculate the SES (standardized effect size) value (Chase, 2010; Zhou et al., 2014). PERMDISP revealed a general pattern in both soil bacterial and archaeal communities, in which ±30% precipitation increased the SES value and ±60% precipitation deceased it (Figures 3C and 3D), and this pattern was robust to increases in the number of rare OTUs included in the analysis (Figure S4). This pattern indicated that ±30% precipitation increased the contribution of deterministic processes relative to stochastic processes, whereas ±60% precipitation increased the contribution of stochastic processes relative to deterministic processes. We further examined the results using a direct-calculation method, and these results were consistent with the changes in the SES values, thus confirming the pattern in Figures 3C and 3D (see details in Data S1, Figure S5, and Table S7).

The increased contribution of deterministic processes relative to stochastic processes under ±30% precipitation was consistent with the traditional view that environmental changes promote deterministic community assembly (Zhang et al., 2011). Both environmental filtering and competitive exclusion are deterministic processes. Although environmental filtering will select for species with similar niches and decrease within-treatment community dissimilarity, interspecific competition can result in the exclusion of species with similar niches and select for species with different niches, thus increasing within-treatment dissimilarity (Chase, 2010; Zhou et al., 2014). Among all five cases where significant differences between the observed within-treatment dissimilarity and the expected dissimilarity of stochastically assembled communities were detected (p < 0.05; Table S6), four cases occurred in the ±30% precipitation treatments, and the observed within-treatment dissimilarity was significantly larger than the expected dissimilarity of stochastically assembled communities (p < 0.05; Table S6), suggesting that ±30% precipitation increased the contribution of the deterministic process of competitive exclusion rather than the environmental filtering process (Chase, 2010; Zhou et al., 2014). Overall, this result supported our hypothesis in Figure 1B that the below-threshold precipitation changes would cause a small increase in the importance of the deterministic process of competitive exclusion in soil microbial community assembly.

The increased contribution of stochastic processes relative to deterministic processes under ±60% precipitation is in contrast to the traditional view that increasing the intensity of environmental changes increases the importance of deterministic processes in community assembly (Zhang et al., 2011). Ecological stochastic processes generally include speciation, species dispersal and colonization, and ecological drift (random birth and death) (Hubbell, 2001). Here, OTUs were defined as 16S rDNA clusters with >97% similarity; although genetic changes might take place in some rapidly evolving genes over the 5-year experimental temporal scale, the 16S rRNA gene was so conservative that evolution could not have had a noticeable effect on the OTU speciation with >97% similarity (Stackebrandt and Goebel, 1994; Zhang et al., 2015). In other words, the increased contribution of stochastic processes is not related to the speciation process. Although the experimental treatments did not affect OTU dispersal, they affected the soil physicochemical conditions (Figure 2) and thus the colonization success of migrating OTUs. The contribution of OTU dispersal/colonization to microbial community assembly was estimated using Hubbell’s neutral model (Hubbell, 2001). The soil bacterial and archaeal OTU migration rate was 0.36–0.40 in these treatments (Table S8), indicating that the stochastic dispersal/colonization process played a role in soil microbial community assembly. However, there were no significant differences in the OTU migration rate among the different precipitation changes (p > 0.05; Table S8), indicating that the stochastic dispersal/colonization process was not the main driver of the increased contribution of stochastic processes. The sharp increase (or decrease) in microbial biomass under +60% (or −60%) precipitation relative to +30% (or −30%) precipitation (Figure 2K) coincided with the increased contribution of stochastic processes under +60% (or −60%) precipitation. In other words, +60% (−60%) precipitation altered the soil habitat in a way that supported the growth of more (or fewer) microbial individuals and resulted in random births (or deaths), which resulted in an increase in the importance of stochastic processes relative to deterministic processes. More advanced technologies and/or refined experimental designs are needed to directly determine the microbial birth/death rates in situ in future studies. Overall, this result supported our hypothesis in Figure 1B that the above-threshold precipitation changes affect soil microbial community assembly by stimulating the stochastic process of ecological drift (random birth/death).

Taxonomic stochasticity translated into functional stochasticity and affected its predictability

To investigate whether changes in the within-treatment taxonomic compositional dissimilarity of soil microbial communities induced by changes in precipitation translated into changes in functional dissimilarities, we characterized the soil microbial metagenome and metaproteome and measured enzyme activities. Consistent with the pattern in taxonomic composition, no systematic responses of the functional composition to precipitation changes were observed (PERMANOVA: metagenome, F = 1.157, p = 0.315; metaproteome, F = 1.038, p = 0.390; enzyme, F = 0.751, p = 0.676); only nonsystematic responses were observed (PERMDISP: metagenome, F = 5.389, p = 0.011; metaproteome, F = 5.988, p = 0.006; enzyme, F = 5.456, p = 0.007; Figures 4A–4C). In particular, the ±30% precipitation treatments resulted in small decreases in the within-treatment dissimilarity of genes, proteins, and enzymes relative to the control, and the ±60% precipitation resulted in large increases in these within-treatment functional dissimilarities (Figures 4A–4C), which were consistent with the changes in within-treatment taxonomic dissimilarity (Figures 3A and 3B). In addition, Mantel tests revealed that the dissimilarities in genes, proteins, and enzymes were significantly correlated with the taxonomic dissimilarities (p < 0.05; Figures 4D–4F) of the soil bacterial community, which accounted for up to 97.64 ± 0.06% of the metagenomic reads in the soil ecosystem (Table S9). Consistent with the significant PERMDISP results (Figures 4A–4C), 24.03% of genes, 25.11% of proteins, and 33.33% of enzymes had relative abundances that were significantly heterogeneous across the five precipitation treatments, and the highest variance in relative abundance often occurred in the +60% or −60% precipitation treatments (p < 0.05; Tables S10, S11, and S12). Furthermore, the within-treatment gene dissimilarity based on the assembled sequences (Figure S6A) and within-treatment dissimilarity based on all the annotated and unannotated reads (Figure S6B; calculated using the Mash [Ondov et al., 2016]) of the soil microbial community showed patterns similar to those of the within-treatment gene dissimilarity based on annotated reads (Figure 4A) and within-treatment taxonomic dissimilarity (Figures 3A and 3B). Overall, these results supported our hypothesis in Figure 1D that the observed change in taxonomic stochasticity translated into functional stochasticity and thus affected functional predictability.

Figure 4

Precipitation changes affect community functions

Effect of precipitation changes on within-treatment dissimilarity of genes, proteins, and enzymes of the soil microbial community revealed by PERMDISP (A–C) and the relationships between genes, proteins, enzymes dissimilarity, and bacteria compositional dissimilarity revealed by the Mantel test (D–F). For multiple comparisons (A–C), the significant (p < 0.05) differences are shown with different uppercase letters. Error bars represent standard errors. See also Figure S6 and Tables S9, S10, S11, and S12.

Structural equation modeling (SEM) was used to determine how soil physicochemical indices, plant community indices, and microbial taxonomic composition were related to microbial community functional indices (Figure S7). The final SEM model adequately fit the data describing the explanatory variable for these functional indices (metagenome, χ2 = 3.24, p = 0.36; metaproteome, χ2 = 3.24, p = 0.52; enzyme, χ2 = 3.24, p = 0.52; standardized path coefficients are given in Figure 5). Among these identified explanatory variables, bacterial taxonomic composition was the variable with the largest absolutes of the standardized path coefficients (0.64–0.94; Figure 5) for all three functional indices. These results suggested that taxonomic information of soil microbial communities is needed to effectively predict the functional responses of soil microbial communities to precipitation changes (Karhu et al., 2014; Zhou et al., 2012).

Figure 5

Explanatory variables for community functions

(A–C) Structural equation model analysis of the effects of soil physicochemical indices, plant community indices, and bacterial taxonomic composition on microbial community function indices (A, metagenome; B, proteome; C, enzyme)

Numbers at solid arrows (p < 0.05) are standardized path coefficients (equivalent to correlation coefficients), and arrow width indicates the strength of the relationships. The dashed arrows indicate nonsignificant relationships (p > 0.05). The number (R2) close to the three functional variables indicate the variance explained by the model. See also Figures S6 and S7 and Tables S10, S11, and S12.

Stochastic process was a distinctive mechanism of above-threshold precipitation changes affecting microbial community assembly

Given that the increased contribution of stochastic processes relative to deterministic processes under ±60% precipitation contrasted with long-standing ideas, we hypothesized that this might represent a novel mechanism by which above-threshold precipitation changes affect soil microbial community assembly. To test this hypothesis, we re-analyzed the soil microbial community composition data of a previous experiment with 13 environmental changes in the same steppe ecosystem (see details of the experimental design in Data S2) (Zhang et al., 2013; Zhou et al., 2020b) to determine whether above-threshold precipitation change had a larger effect on bacterial community size and stochastic ecological processes compared with other environmental changes. In the increased precipitation treatment in the multifactorial experiment, the precipitation was approximately +54.43% higher compared with natural levels of precipitation during the plant growing season (from May to August) over the 5-year treatment. Changes in each environmental factor might have unique effects on the soil microbial community due to their distinct properties; for example, climate warming and nutrient enrichment affect microbial communities by increasing the soil temperature and nutrient content, respectively. Environmental changes can also vary in intensity, and variation in the intensity of environmental changes can affect the magnitude of deterministic changes in community composition. To control the effects of differences in the intensity of the 13 environmental changes and characterize the effects of their distinct properties, we standardized the stochastic change in bacterial community composition (which was calculated as described in Data S1) using the following formula: the absolute value of ([the stochastic change in community composition]/[the deterministic change in community composition]). The change in soil bacterial community size was standardized using the following formula: the absolute value of ([percentage change in community size relative to the control]/[the deterministic change in community composition]) (Zhang et al., 2011, 2016). The standardized change in bacterial community size was positively correlated with the standardized stochastic change in bacterial community composition (slope = 0.056, F = 35.420, p < 0.001, R2 = 0.763; Figure 6B), and this positive linear relationship remained even when the two points with the largest values of standardized change in bacterial community size (watering, watering + warming) were excluded from the analysis (slope = 0.032, F = 4.658, p = 0.059, R2 = 0.341). These results indicated that the above-threshold precipitation increase caused the soil habitats to support more microbial individuals (larger community size) and thus stimulated the random birth process, which increased the relative contribution of stochastic processes to microbial community assembly. Among all the five precipitation-associated environmental changes, there were three with relatively large values of standardized stochastic change in bacterial community composition (watering, watering + warming, and mowing + watering; Figure 6A), demonstrating that the above-threshold precipitation increase had a greater effect on bacterial community size and the stochastic birth process compared with other factors. The other two precipitation-associated environmental changes with relatively small values of standardized stochastic change in bacterial community composition were both N-associated (N addition + watering, mowing + N addition + watering; Figure 6A), suggesting that the stimulating effect of precipitation increase was counteracted by N addition. N addition has often been observed to decrease soil bacterial community size and alter bacterial community composition by reducing soil pH in this steppe ecosystem (Zhang et al., 2011, 2015); thus, N addition and watering had opposite effects on bacterial community size, which indicates that the effects of these treatments were counteracted when they were combined, such as in the “N addition + watering” treatment. The standardized stochastic change in bacterial community composition was also high for warming (Figure 6A), the ecological consequence of which requires further study. Anyway, these results showed that the above-threshold precipitation increase and its co-occurrence with other factors affected soil microbial community assembly through the distinctive process of stochastic birth.

Figure 6

Evidence for the distinctive pattern and mechanism of precipitation changes affecting soil microbial assemblages

(A and B) Standardized changes in soil bacterial community size under 13 anthropogenic environmental changes (A) and their relationship with the standardized stochastic changes in bacterial community composition (B). See also Data S2.

Conclusions

We revealed the distinctive patterns and mechanisms by which precipitation changes affect soil microbial community assembly in Eurasian steppe and confirm all components of our hypothesis (Figure 1). Briefly, below-threshold precipitation increases (or decreases) increase the importance of the deterministic process of competitive exclusion in soil microbial community assembly and thus decrease the within-treatment compositional dissimilarity, which results in decreased functional dissimilarity. This will cause the community functions to become more predictable with the soil physicochemical indices alone, making the microbial taxonomic composition data unnecessary for the prediction. In contrast, above-threshold precipitation increases (or decreases) stimulate the random birth (or death) process; this leads to increases in within-treatment compositional dissimilarity and causes community functions to become more unpredictable, which necessitates the use of taxonomic composition data.

Our focus in this study was on community functional composition, such as the relative abundances of different SOM (soil organic matter) catabolic pathways (Figure S8); thus, the taxonomic stochasticity led to functional stochasticity and increased its unpredictability. If we were focused on overall community functions (e.g., community respiration) or individual functions (e.g., N fixation), the size of the functions would be more predictable. For example, the amount of community respiration might show a linear trend with the precipitation gradient, which is similar to the pattern of microbial community size (Figure 2K). The effects and mechanisms of precipitation changes on various overall functions and special functions require further study.

How natural precipitation gradients across different ecosystems (e.g., from the equator to the Arctic) affect the relative contributions of stochastic processes and deterministic processes to soil microbial community assembly also requires additional study (Fierer, 2017). Here, we hypothesized that the random birth process would lead to more mutation and even speciation, and the random death process would lead to species extinction over large spatiotemporal scales. The stochastic effects and phenomena are more vulnerable to being neglected; thus, the distinctive stochastic mechanism of high-intensity precipitation change suggested that the importance of precipitation in driving soil microbial assemblages has likely been underestimated in previous studies. In particular, the effect of precipitation or soil water content on microbial community assembly may be partly covered by soil pH. For example, in the control communities of this study, there was a significant positive and negative relationship of archaeal community composition PCoA (principal coordinate analysis) axis 1 with the soil water content and pH, respectively (Table S7), which indicated that soil environmental heterogeneity at a small spatial scale caused variation in microbial community composition. Thus, teasing apart the effect of soil pH is important for characterizing the stochastic mechanism by which precipitation changes affect soil microbial assembly.

Limitations of the study

In this study, ±30% and ±60% precipitation were used as the below- and above-threshold precipitation changes, respectively. Experiments with multiple precipitation treatments are needed to determine the exact thresholds where deterministic processes give way to stochastic processes, and the thresholds might be different for precipitation increases and decreases (Zhang et al., 2017a). The thresholds likely have special physiological implications; here, we hypothesized that microorganisms can adapt to consuming added resources (or tolerate resource shortages) in response to below-threshold precipitation increases (or decreases) and that they are unable to adapt physiologically in response to above-threshold precipitation increases (or decreases), but they can consume added resources by increasing their abundances through the birth process (or deal with a shortage of resources by decreasing their abundances; i.e., the death process). In addition, more time-series samples are needed to monitor the random birth/death processes in future studies.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ximei Zhang (zhangximei@caas.cn).

Materials availability

This study did not generate new unique reagents.

Data and code availability

All amplicon and metagenome sequences have been deposited in NCBI, and metaproteome sequences have been deposited in the ProteomeXchange Consortium, with all these accession numbers listed in the method section. Any additional information required to re-analyze the data reported in this paper are available from the lead contact upon request.

Method details

Study site and experimental design

This study was conducted in a semiarid steppe ecosystem (44°22ʹ N, 117°35ʹ E) located in West Ujinmqin Banner, Inner Mongolia, China. The site has not experienced grazing or other anthropogenic disturbances since it was fenced in 2011. Before being fenced, the ecosystem was mown once a year in August. The mean annual air temperature was ∼1.6°C, and the annual precipitation was ∼334 mm, with ∼86% of the precipitation falling during the growing season (May–September). The site is characterized by typical steppe vegetation dominated by C3 perennial grasses and forbs, such as Anemarrhena asphodeloides, Leymus chinensis, and Stipa grandis.

The experiment was conducted in a randomized block design with five precipitation levels: −0.6 (precipitation decreased by 60%), −0.3 (precipitation decreased by 30%), 0 (ambient control; natural precipitation), 0.3 (precipitation increased by 30%), and 0.6 (precipitation increased by 60%). The historical growing-season precipitation varied from 59% below to 77% above the long-term average precipitation in the area; thus, the ±0.6 and ±0.3 levels were used to represent the above- and below-threshold precipitation changes, respectively. A total of 20 plots were assigned to four blocks, and each block contained the five randomly assigned treatment plots. Each plot was 3 × 4 m with 1 m of spacing between plots. To limit belowground water movement between plots, we inserted tin sheets into the ground to a 100 cm depth around each plot. There was a 0.5 m buffer on each side of the plot, and all measurements were conducted in the central area (2 × 3 m) to avoid edge effects. Precipitation manipulations were carried out from May (June in 2012) to August during 2012–2017. In the −0.6 and −0.3 plots, 60% and 30% of the rainfall, respectively, was blocked using passive rainout shelters according to a previously published methodology (Yahdjian and Sala, 2002). Although the shelters intercept incoming light to a small degree, we have previously shown that this interception of light has a negligible influence on air temperature and plant photosynthesis (Zhang et al., 2017a). In the 0.3 and 0.6 plots, 30% and 60% rainfall, respectively, was added immediately to the plots after each rainfall event (>2 mm) by a handheld irrigation system with a flow meter. All added water was rainwater collected from the shelters.

Sampling and measurement of soil physicochemical indexes

The volumetric water content and temperature of surface soil (0–10 cm) in each plot were monitored every 10 days during the growing season with a TDR-300 soil moisture probe (Spectrum Technologies Inc., Plainfield, Illinois, USA) and a thermocouple probe connected to an LI-8100 (LI-COR Inc., Lincoln, Nebraska, USA). In mid-August 2017 (the period when plant biomass is highest), aboveground vegetation was sampled by clipping all plants at the soil surface from a 0.2 × 0.8 m strip randomly placed in the plot. All living vascular plants were sorted into species, oven-dried at 65ºC for 48 h, and weighed. The number of species detected was used to represent plant richness, and the dry mass of all living plants was used to represent the aboveground plant biomass and community composition.

In mid-August 2017, four soil cores (10 cm deep, 3.5 cm diameter) were collected from each plot at random and thoroughly homogenized. After removing roots and stones with a 2 mm sieve, part of the soil sample was frozen in liquid N for measurement of microbial meta-proteins; another part was frozen at −20°C for DNA extraction; another part was kept fresh at 4°C for measurement of the dissolved organic C content, inorganic available N (NH4+-N and NO3--N) content, microbial biomass C and N content, soil enzyme activities; and the remaining portion was kept dry for measurement of soil pH, organic C content, and total N content.

Dissolved organic C was extracted by adding 50 ml of 0.5 M potassium sulfate (K2SO4) to subsamples of 12.5 g of homogenized soil and agitating them on an orbital shaker at 120 rpm for 1 h. The filtrate was analyzed using a TOC analyzer (HighTOC, Elementar). Soil NH4+-N and NO3--N content was determined on a FIAstar 5000 Analyzer (Foss Tecator, Denmark) after extraction of fresh soil with 1 mol/L KCl within a month of soil collection (Bao, 2000). Soil pH was measured in 1:2.5 (W/V) suspensions of soil in distilled water. The soil organic C content and total N content were quantified using the potassium dichromate-vitriol oxidization method and the Kjeldahl acid-digestion method, respectively (Bao, 2000).

Soil microbial biomass C and N content was estimated using the chloroform fumigation–extraction method (Vance et al., 1987). Briefly, aliquots of the fresh soil (15 g dry weight equivalent) were fumigated for 24 h with ethanol-free CHCl3. Additional aliquots of fresh soil were used as unfumigated controls. Both the fumigated and unfumigated samples were extracted with 60 mL of 0.5 M K2SO4 for 30 min on a shaker. K2SO4 extracts were filtered through 0.45 mm filters and frozen at −20°C before analysis of extractable C and N with a Liqui TOC elementar analyzer (Elementar Liqui TOC, Elementar Co., Hanau, Germany). The microbial biomass C and N content was calculated from differences in the extractable C and N content between the fumigated and unfumigated samples using conversion factors (kEC and kEN) of 0.45.

Measurement of soil enzyme activities

Soil enzyme activities were measured following the procedures of DeForest (2009). Buffer and substrate recipes followed those of Saiya-Cork et al. (2002) but with several modifications. A 0.1-mM acetate buffer solution was made by mixing sodium acetate trihydrate (certified ACS grade, crystalline, Fisher Scientific) with deionized (DI) water. The pH was adjusted using 12 M HCl (technical grade, Fisher Scientific) and stored in a 20−1 HDPE carboy. Because enzyme activity is pH sensitive, the pH of the buffer was adjusted to within 0.5 units of the mean soil pH of the study site. Acetate buffer was kept at 4°C between analyses and used within 8 days; new buffer was made on days 14 and 21.

Activities of the β-glucosidase, β-xylosidase, acid phosphatase, cellobiohydrolase, and NAGase were measured fluorometrically using MUB-linked model substrates, and the activity of leucine aminopeptidase was measured fluorometrically using AMC-linked model substrate (Marx et al., 2001; Saiya-Cork et al., 2002) (Table S13). Each MUB-linked substrate (or AMC-linked substrate) (Sigma–Aldrich, St. Louis, MO) was a 100-mM solution mixed well in DI water and stored in opaque containers.

Approximately 1.00 ± 0.01 g of soil was placed in a 400 ml amber HDPE bottle; time was recorded, and samples were kept cool until all the samples were weighed. The time between when the soil was weighed and buffer was added was approximately 22 min. Approximately 100 ml of the 0.1 mM sodium acetate buffer was added, and the samples were homogenized with a magnetic stirrer at high speed for 1 min.

Activities of these enzymes were measured fluorometrically in black polystyrene 96-well, 300-μl microplates (Whatman Inc., Florham Park, NJ). Buffer, samples, references, and substrates were dispensed in a strict order and position on the well plate. First, 200 μl of buffer was dispensed in columns that served as the blank, reference, and negative control. Next, 50 μl of buffer was dispensed in columns that served as the blank and sample controls for three soils. The purpose of the blank and sample controls was to account for the naturally occurring fluorescence of the buffer or soil within the well plate. Next, 200 μl of the blended soil slurry was dispensed into columns that served as the quench, soil control, and sample assay. Fifty μl of a 100-μM MUB (or AMC-linked substrate) solution was dispensed into columns that served as a reference standard and quench standard. Lastly, 50 μl of a 100-μM MUB-linked substrate (or AMC-linked substrate) was dispensed into columns that served as the negative control and sample assays. The addition of the MUB-linked substrate (or AMC-linked substrate) marked the start of the incubation period, and the time was noted. The microplates were covered and incubated in the dark at 20°C for 0.5 h for N-acetyl-glucosaminidase and phosphatase and 2 h for β-glucosidase, β-xylosidase, and cellobiohydrolase (Saiya-Cork et al., 2002). The fluorescence was measured using a microplate fluorometer (Synergy HT, BioTek, Winooski, VT) with 365-nm excitation and 450-nm emission filters (Saiya-Cork et al., 2002). Enzyme activities were expressed in units of nmol h−1 g−1 and calculated using the following equations:where

The 100 μl refers to the volume of soil slurry, and 0.5 nmol is the amount of MUB standard (or AMC standard) added to a well. Fluorescence values are means from the eight analytical wells subtracted from the blank. Soil (g) is the oven-dried weight.

Soil microbial taxonomic diversity measurement

Microbial community genomic DNA was extracted from soil samples using the MoBio PowerLyzer PowerSoil DNA isolation kit according to manufacturer’s instructions. The DNA extract was checked on 1% agarose gel, and DNA concentration and purity were determined with NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA).

The hypervariable region V3-V4 of the bacterial 16S rRNA gene was amplified using an ABI GeneAmp® 9700 PCR thermocycler (ABI, CA, USA) with the primer pairs 338F (5′-ACT CCT ACG GGA GGC AGC AG-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′), and the hypervariable region V4 of the archaeal 16S rRNA gene was amplified with the primer pairs 524F10extF (5′-TGY CAG CCG CCG CGG TAA-3′) and Arch958RmodR (5′-YCC GGC GTT GAV TCC AAT T-3′). The PCR mixtures contained 4 μL of 5 × TransStart FastPfu buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of forward primer (5 μM), 0.8 μL of reverse primer (5 μM), 0.4 μL of TransStart FastPfu DNA Polymerase, 10 ng of template DNA, and up to 20 μL of ddH2O. The PCR amplification was performed as follows: initial denaturation at 95°C for 3 min; 28 (for bacteria) or 36 (for archaea) cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 45 s; and a single extension at 72°C for 10 min, ending at 4°C. PCR reactions were performed in triplicate. The PCR products were extracted from a 2% agarose gel, purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instructions, and quantified using a Quantus™ Fluorometer (Promega, USA).

Purified amplicons were pooled in equimolar amounts and paired-end sequenced on an Illumina MiSeq PE300 platform (Illumina, San Diego, USA) according to the standard protocols of Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: SRP295607 for bacteria, SRP295606 for archaea).

The raw 16S rRNA gene sequencing reads were demultiplexed, quality-filtered by fastp version 0.20.0 (Chen et al., 2018), and merged by FLASH version 1.2.7 (Magoc and Salzberg, 2011) using the following criteria: (i) 300 bp reads were truncated at any site receiving an average quality score of <20 over a 50 bp sliding window, truncated reads shorter than 50 bp were discarded, and reads containing ambiguous characters were discarded (Table S14); (ii) only overlapping sequences longer than 10 bp were assembled, the maximum mismatch ratio of the overlap region was 0.2, and reads that could not be assembled were discarded; and (iii) samples were distinguished according to the barcode (exact barcode matching, barcode length was 10 bp, and no mismatches were allowed) and primers (two nucleotide mismatches were permitted for primer matching), and the sequence direction was adjusted (Table S14).

OTUs (97% similarity cutoff) were clustered using UPARSE version 7.1 (Edgar, 2013; Stackebrandt and Goebel, 1994), and chimeric sequences were identified and removed. In addition, singleton reads were discarded. The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 (Wang et al., 2007) against the 16S rRNA database (eg. Silva v138) using a confidence threshold of 0.7.

Metagenomic sequencing and analysis

The extracted DNA was fragmented to an average size of approximately 260 bp using Covaris M220 (Gene Company Limited, China) for paired-end library construction. The paired-end library was constructed using NEXTFLEX Rapid DNA-Seq (Bioo Scientific, Austin, TX, USA). Adapters containing the full complement of sequencing primer hybridization sites were ligated to the blunt-end of fragments. Paired-end sequencing (paired-end 150 bp) was performed on an Illumina Hiseq Xten platform (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using HiSeq X Reagent Kits according to the manufacturer’s instructions (www.illumina.com). Sequence data have been deposited in the NCBI Short Read Archive database (Accession Number: SRP295741). Shotgun sequencing resulted in 80.15 ± 1.13 million sequences (mean ± one standard error) or 11.27 ± 0.16 Giga base pairs (Gbp) of sequences per sample.

Paired reads of the shotgun metagenomic sequences were merged using FLASH with default parameters (Magoc and Salzberg, 2011). To improve the reliability and quality of the subsequent analysis, the raw sequence data were processed using Seqprep software (https://github.com/jstjohn/SeqPrep) to remove the adapter sequences. Next, the library sickle (https://github.com/najoshi/sickle) was used to trim the reads from the 5′ end to the 3′ end using a sliding window (size 50 bp, 1 bp steps) (Joshi and Fass, 2011). Sequences were trimmed where the mean quality of bases inside a window was below 20 (99% accuracy), and the longest contiguous sequence above this threshold was retained for further analysis. Quality-trimmed reads shorter than 50 bp or containing N (ambiguous bases) were discarded.

Metagenomics data were assembled using MEGAHIT (Li et al., 2015) (https://github.com/voutcn/megahit, version 1.1.2), which makes use of succinct de Bruijn graphs. The final assembly results consisted of contigs with lengths 300 bp or greater.

Open reading frames (ORFs) from each merged read (or assembled contig) were predicted using MetaGene (http://metagene.cb.k.u-tokyo.ac.jp/) (Noguchi et al., 2006). The predicted ORFs with lengths 100 bp or greater were retrieved and translated into amino acid sequences using the NCBI translation table (http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes#SG1).

All predicted genes with a 95% sequence identity (90% coverage) were clustered using CD-HIT (http://www.bioinformatics.org/cd-hit/) (Fu et al., 2012), the longest sequences from each cluster were selected as representative sequences to construct a non-redundant gene catalog. Reads after quality control were mapped to the representative sequences with 95% identity using SOAPaligner (http://soap.genomics.org.cn/) (Li et al., 2008), and the gene abundance in each sample was evaluated.

Representative sequences of the non-redundant gene catalog were aligned to the NCBI NR database with an e-value cutoff of 1e-5 using BLASTP (Version 2.2.28+, http://blast.ncbi.nlm.nih.gov/Blast.cgi) for taxonomic annotations (Altschul et al., 1997). Cluster of orthologous groups of proteins (COG) annotations for the representative sequences were determined using BLASTP against the eggNOG database (Version4.5) with an e-value cutoff of 1e-5 (Jensen et al., 2008; Tatusov et al., 2003).

The relative abundance of each COG gene was determined by the sum of the reads mapped to each COG gene (Tatusov et al., 2000) normalized by the size of the dataset. We specifically focused on the effect of precipitation changes on the abundances of genes for the catabolic processes of various C complexes differing in decomposability, ranging from the highly recalcitrant lipids and phenolics to the more labile monosaccharides, sugar acids, and sugar alcohols (Zhang et al., 2017b).

Mash, a tool that uses kmers to compare the sequence composition between metagenomes, was used to calculate sample-sample dissimilarities (options: -k 25 –s 100,000; Ondov et al., 2016).

Meta-protein measurements

Samples were ground with liquid N, and BPP was added in a ratio of 1:10. The solution was centrifuged at 12,000 g for 20 min at 4°C, and the supernatant was collected. An equal volume of Tris-saturated phenol was added and vortexed for 10 min at 4°C. The solution was centrifuged at 12,000 g for 20 min at 4°C, and the phenol phase was collected. An equal volume of BPP was added and vortexed for 10 min at 4°C. The solution was centrifuged at 12,000 g for 20 min at 4°C, and the phenol phase was collected. Five volumes of pre-cooled 0.1 M ammonium acetate in methanol were added, and protein was precipitated overnight at −20°C. The supernatant was discarded by centrifugation, and the precipitate was washed twice with 90% acetone. After discarding the supernatant by centrifugation, the precipitate was re-suspended with lysis buffer (1% SDS, 8 M urea, cocktail) and then sonicated for 2 min on ice. Finally, the lysates were centrifuged at 12,000 g for 20 min at 4°C, and the supernatants were collected to determine the concentration of protein in all samples. Protein concentrations were determined by the bicinchoninic acid (BCA) method by a BCA Protein Assay Kit (Thermo Fisher Scientific). Protein quantification was performed according to the kit protocol.

Protein digestion was performed according to the standard procedure. For each sample tube containing 200 μg of protein, TCEP was added to a final concentration of 10 mM, and the tubes were incubated at 37°C for 60 min. IAM was then added to a final concentration of 40 mM and was left to react for 40 min in the dark. Five volumes of cold acetone were added to the sample tube. The tube was then inverted three times and incubated at −20°C until the precipitate formed (∼4 h). The acetone was removed by centrifugation at 10,000 g for 20 min, and the precipitated protein was resuspended with 100 mM TEAB Buffer. Trypsin solution was added to each sample tube at a proportion of 1:50, and the tubes were incubated at 37°C overnight.

The peptides were vacuum-dried and then resuspended with 2% acetonitrile and 0.1% TFA. Samples were desalted with HLB and MCX and vacuum-dried. Peptide concentrations were determined by a peptide quantification kit (Thermo, Cat.23275). Loading buffer was added to each tube to prepare samples for mass spectrometry analysis, and the concentration of each sample was 0.5 μg/μl.

Experiments were performed on a Q-Exactive that was coupled to an Easy-nLC 1200 system. Each peptide sample was injected for nanoLC-MS/MS analysis. The sample was loaded onto a C18 reversed-phase column (75 μm × 25 cm, Thermo Fisher Scientific, USA) in buffer A (2% acetonitrile and 0.1% formic acid) and separated with a linear gradient of buffer B (80% acetonitrile and 0.1% formic acid). The peptides were eluted using the following gradient: 0–2 min, 0–5%B; 2–105 min, 5–23%B; 105–130 min, 23–29%B; 130–147 min, 29–38%B; 147–149 min, 38–100%B; 149–155 min, 100−100%B; 155–160 min, 100−0%B; and 160–180 min, 0−0%B; the flow rate was 300 nL/min.

An electrospray voltage of 1.8 kV versus the inlet of the mass spectrometer was used. Q-Exactive was operated in the data-dependent mode to switch automatically between MS and MS/MS acquisition. Survey full-scan MS spectra (m/z 350–1300) were acquired with a mass resolution of 70K, followed by twenty sequential high-energy collisional dissociation (HCD) MS/MS scans with a resolution of 17.5K. In all cases, one microscan was recorded using a dynamic exclusion of 30 s.

MS/MS spectra were searched using ProteomeDiscovererTM Software 2.2 against an already established database. The highest score for a given peptide mass (best match to that predicted in the database) was used to identify parent proteins. The parameters for protein searching were set as follows: tryptic digestion with up to two missed cleavages, carbamidomethylation of cysteines as fixed modification, and oxidation of methionines and protein N-terminal acetylation as variable modifications. Peptide spectral matches were validated based on q-values at a 1% false discovery rate (FDR).

A total of 4,062 proteins belonging to the proteome of the soil samples were identified in this study. All identified proteins were annotated using Gene Ontology (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes pathway (http://www.genome.jp/kegg/) analysis. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the iProX partner repository with the dataset identifier PXD023103.

Relative contribution of deterministic and stochastic processes

β diversity, which represents the compositional variation between communities, is often used to infer mechanisms of community assembly, such as the relative importance of deterministic versus stochastic processes. However, differences in β-diversity indexes may be caused by differences in ecological processes, as well as α and γ diversity. To account for the effect of the other two diversity components, Chase (2010) developed a null model method that compares the observed β diversity to the theoretical β diversity from stochastically assembled communities. Because this method depends on the presence or absence of OTUs and is sensitive to noise from rare OTUs, various relative abundance thresholds for OTUs (>1%, 0.1%, 0.01%, 0.001%, 1/148,106, or 1/716,472; Table S1) were used in analyses (Ferrenberg et al., 2013). To identify the relative contributions of deterministic and stochastic processes to driving soil microbial assembly, we analyzed the community data of each treatment following the methods of Chase (2010) and Zhou et al. (2014). First, for any given pair of plots within a treatment, we calculated the observed OTU richness (e.g., α1 and α2 for plot 1 and 2, respectively) and the number of shared OTUs (SSobs). Second, the total number of OTUs detected in the “OTU pool” (γ diversity) from all plots of a treatment and the proportion of the plots occupied by each OTU were measured. Third, we calculated the distribution of the expected shared OTUs from the null model (SSexp) by randomly drawing α1 and α2 OTUs from the OTU pool, and the probability of an OTU drawn was proportional to its among-plot occupancy. The SSexp and the expected Jaccard’|'s similarity (Jexp) were obtained for each draw, and the average Jaccard’s similarity () and its SD were estimated based on 10,000 draws (σexp). For each treatment, PERMDISP was used to test the difference between the observed community similarity (Jobs) and the average of the expected community similarity (). Non-significant differences (p > 0.05) indicated that stochastic processes were the primary driver of community assembly and significantly larger (or smaller) Jobs relative to indicates that the deterministic process of ecological filtering (or competitive conclusion) was the primary driver of community assembly.

To investigate whether the treatments affected the relative contributions of deterministic and stochastic processes to driving microbial community assembly, the community data of all plots of all treatments were used to calculate Jobs and (and its σexp) following the steps described above. Thus, the OTUs of all treatments represented the OTU pool. The magnitude of the effect of deterministic processes on community structure was further quantified using the SES index (Chase, 2010; Zhou et al., 2014): SES = (Jobs–)/σexp. PERMDISP was used to reveal the effect of the precipitation gradient on the SES value.

We also used a direct-calculation method to parse the relative importance of the deterministic and stochastic components of the treatment effect following previously described methods (Zhang et al., 2011, 2016). Formulas for the deterministic change and stochastic change of each treatment are as follows: deterministic change of each treatment = ((mean compositional variation between treatment and control) – (mean compositional variation within control)); stochastic change of each treatment = ((mean compositional variation within treatment) – (mean compositional variation within control)).

Quantification and statistical analysis

Most statistical analyses were conducted in R unless otherwise stated. To account for the effect of unequal sampling in the amplicon sequencing and metagenomic sequencing on the various community indices, bacterial OTUs with relative abundance <1/716472 and archaeal OTUs with relative abundance < 1/148106 were removed in the analyses (716,472 and 148,106 are the smallest sequence numbers for bacteria and archaea, respectively). Bray–Curtis dissimilarity based on the relative abundance of OTUs was calculated and used to characterize variation in community composition among samples using the function “vegdist” in the “vegan” package (Oksanen et al., 2020), and PCoA was used to visualize the differences among these samples with the function “pcoa” in the package “ape” (Paradis and Schliep, 2019). The Bray–Curtis dissimilarity based on relative abundances of genes and proteins was also calculated. For the soil enzyme activity data, we first standardized the data for different samples from zero to one and then calculated the Bray–Curtis dissimilarity based on the data of all the enzymes. PERMANOVA and PERMDISP were used to reveal the effect of the precipitation gradient on the Bray–Curtis dissimilarity matrices with the functions “adonis” and “betadisper” in the package “vegan” (Oksanen et al., 2020), respectively. Pearson correlations were used to determine the relationships between the relative abundances of each OTU (or gene, protein, and enzyme) and the soil water content, and the p-values were corrected using the FDR method and the function “corr.test” with “adjust = fdr” in the “psych” package (Revelle, 2020). Levene’s test was used to test whether the variances in the relative abundance of each OTU (or gene, protein, and enzyme) were homogeneous among the five precipitation gradients with “leveneTest” in the package “car” (Fox and Weisberg, 2019). A Mantel test was used to assess the relationships between microbial community composition and the plant and soil indices, and the relationships between the Bray-Curtis dissimilarity matrices were analyzed using the function “mantel” in the package “vegan” (Oksanen et al., 2020).

The variation partitioning analysis (VPA) was performed to quantify the relative contributions of soil factors (soil temperature, pH, water content, organic C content, total N content, dissolved organic C content, and available N content), plant biomass, and plant community composition to bacteria and archaea community composition. We conducted a co-linearity test to identify co-linear factors; the soil organic C content and total N content were removed from subsequent analyses because their variance inflation factor values were >10. The first two axes of the PCoA were used to characterize plant community composition. The dissimilarity-based redundancy analysis (db-RDA) and forward selection procedure based on the Bray-Curtis dissimilarity matrices were used to select the subsets of soil factors (soil temperature, pH, water content, dissolved organic C content, and available N content) using the "dbrda" and "ordistep" functions in the "vegan" package (Oksanen et al., 2020). The forward selection was stopped if the p-value was greater than 0.1 or if the Akaike information criterion (AIC) did not decrease after more variables were added. VPA was performed with an adjusted R2 value to determine the effects of selected soil factors, plant biomass, and plant community composition on the bacterial and archaeal communities using the function "varpart" in the "vegan" package (Oksanen et al., 2020). The statistical significance of each group of explanatory variables was evaluated using a permutation test with the function "anova" in the package “stats” (R Core Team, 2020).

Hubbell’s neutral theory for the maintenance of biodiversity in a local community (Hubbell, 2001) predicts that the species abundance distribution at the local scale is a zero-sum multinomial distribution (ZSM) and that the migration rate (m) is the key parameter (Volkov et al., 2003). The parameter m defines the probability that a newly established individual is an immigrant from the outside (i.e., the regional species pool). Dispersal limitation occurs when m is small. Therefore, the migration rate indicates the dispersal limitation by fitting a ZSM by maximum likelihood using the OTU abundance data for each sample, which was defined as the local community (McGill et al., 2006). Specifically, TeTame based on the sampling formula developed by Etienne (2005) has been used to calculate the migration rate (Jabot et al., 2008).

We used structural equation modeling (SEM) to determine how soil physicochemical indices, plant community indices, and bacterial taxonomic composition affected microbial community functional indices (Figure S7). SEM is based on a simultaneous solution procedure, in which the residual effects of predictors are estimated (partial regressions) once common causes from inter-correlations have been statistically controlled for (Grace et al., 2012). The first PCoA axes were used in the SEM analysis to indicate plant community composition, microbial taxonomic and functional composition, and soil variable dissimilarity (soil temperature, pH, dissolved organic C content, and available N content) based on the Bray-Curtis distance. In the SEM analysis, we compared the model-implied variance–covariance matrix against the observed variance–covariance matrix, and the data were fitted to the models using the maximum likelihood estimation method. The adequacy of the models was determined using χ2 tests, and adequate model fits were indicated by a non-significant χ2 test (p > 0.05) (Grace et al., 2012). SEM was conducted by establishing the a priori models shown in Figure S7. Non-significant paths with the largest p values were eliminated in a new model until all the pathways were significant. The best-fit model was selected based on the chi-square test, root mean square error of approximation (RMSEA) value, and AIC. The p values of the chi-square test, RMSEA, and AIC of all sub-models are shown in Table S15. The SEM-related analyses were performed using the lavaan package in R 4.0.2 (R Core Team, 2020; Rosseel, 2012).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Critical commercial assays
MoBio PowerLyzer PowerSoil DNA isolation kit	Mo Bio	Cat. 12855
Deposited data
Amplicon data of 16S rDNA of soil bacterial and archaeal communities	This paper	NCBI Sequence Read Archive database: SRP295607 and SRP295606
Metagenome data of soil microbial community	This paper	NCBI Short Read Archive database: SRP295741
Metaproteome data of soil maicrobial community	This paper	ProteomeXchange Consortium: PXD023103
Oligonucleotides
338F	Majorbio Bio-Pharm Technology Co. Ltd.	5'-ACT CCT ACG GGA GGC AGC AG-3'
806R	Majorbio Bio-Pharm Technology Co. Ltd.	5'-GGA CTA CHV GGG TWT CTA AT-3'
524F10extF	Majorbio Bio-Pharm Technology Co. Ltd.	5'-TGY CAG CCG CCG CGG TAA-3'
Arch958RmodR	Majorbio Bio-Pharm Technology Co. Ltd.	5'-YCC GGC GTT GAV TCC AAT T-3'
Software and algorithms
fastp version 0.20.0	Chen et al., 2018	https://github.com/OpenGene/fastp
FLASH version 1.2.7	Magoc and Salzberg, 2011	http://ccb.jhu.edu/software/FLASH/
UPARSE version 7.1	Edgar, 2013	http://drive5.com/uparse
RDP Classifier version 2.2	Wang et al., 2007	http://rdp.cme.msu.edu/
the library sickle	Joshi and Fass, 2011	https://github.com/najoshi/sickle
MEGAHIT version 1.1.2	Li et al., 2015	https://github.com/voutcn/megahit
MetaGene	Noguchi et al., 2006	http://metagene.cb.k.u-tokyo.ac.jp/
CD-HIT	Fu et al., 2012	http://www.bioinformatics.org/cd-hit/
SOAPaligner	Li et al., 2008	http://soap.genomics.org.cn/
BLASTP version 2.2.28+	Altschul et al., 1997	http://blast.ncbi.nlm.nih.gov/Blast.cgi
Mash	Ondov et al., 2016	https://github.com/marbl/mash
Vegan (R package)	Oksanen et al., 2020	=https://CRAN.R-project.org/package=vegan
ape (R package)	Paradis and Schliep, 2019	http://cran.r-project.org/package=ape
psych (R package)	Revelle, 2020	https://CRAN.R-project.org/package=psych
car (R package)	Fox and Weisberg, 2019	https://CRAN.R-project.org/package=car
stats (R package)	R Core Team, 2020	https://www.R-project.org/
Other
Illumina sequencing	Illumina	MiSeq; Hiseq Xten