Warming and N amendment experiments were carried out at the Sanjiangyuan Alpine Grassland Ecosystem Field Observation Station (N33°24′30″, E97°18′00″, altitude 4,270 m) on Southern Tibetan Plateau, Qinghai, China. This area is an alpine grassland, dominated by Kobresia pygmaea, Lagotis ramalana, Potentilla nivea, and Thalictrum alpinum (Wen and Lu, 2015; Wang et al., 2019). As described previously (Wang et al., 2019), this site is characterized by a continental monsoon climate, comprised of cold (October to May) and growing (June to September) seasons, with an annual mean air temperature from -5.6 to 3.8°C during 1985–2004, and annual average precipitation of 562.2 mm, with ∼75% of annual precipitation occurring during the growing season.

Experimental plots were set up in July 2013, in which terrain was flat and plant community composition was relatively uniform. It was surrounded by an enclosure to avoid the grazing disturbance. Warming experiments were carried out using 4 m² × 0.4 m tall Open-top chambers (OTCs) (Henry and Molau, 1997), constructed with clear polyvinyl chloride (PVC) and a steel framework (Wang et al., 2019). OTCs increased mean near-surface temperatures by 2.00 ± 0.24°C. Meanwhile, either ammonium sulfate or sodium nitrate was amended as the external N source for alpine grassland at 300 kg N ha year^–1 for maximal plant biomass (Guo et al., 2014). Soil samples were collected from the 0–15 to 15–30 cm strata, due to the majority of grass root growing between 0 and 15 cm. As our focus was on the effects of warming on different soil strata, we set group treatments, including ammonium sulfate and sodium nitrate amendment, with nine replicates per group treatment as: (i) warming upper (WUp), warming treatment at 0–15 cm strata; (ii) warming lower (WLo), warming treatment at 15–30 cm strata; (iii) unwarming upper (UWUp), ambient temperature at 0–15 cm strata; and (iv) unwarming lower (UWLo), ambient temperature at 15–30 cm strata.

Soil samples were collected in August 2016. Soil cores, 5 cm in diameter, were taken from 0 to 30 cm at three randomly chosen locations within each 4 m² plot and divided into 0–15 and 15–30 cm layers, then mixed together, respectively, to form composite samples. Soil samples were set on ice then transported back to the laboratory where visible grassroots and pebbles were subsequently removed. Soil samples were stored at −80 and 4°C for DNA extraction and soil geochemical measurements, respectively. Soil geochemical measurements, including pH, ammonium and nitrate content, total nitrogen content, organic matter, available nitrogen content, soil temperature, and soil moisture were described previously (Wang et al., 2019).

Genomic DNA was extracted from 2 g of soil with the freeze-grinding method as described previously (Zhou et al., 1996), and purified using the PowerSoil ^® DNA Isolation Kit (MoBio Laboratories, Inc., Solana Beach, CA, United States) according to the manufacturer’s instructions. For polymerase chain reaction (PCR) amplification, DNA quality and concentration were assessed based on absorbance ratios of 260/280 nm (∼1.8) and 260/230 nm (>1.7) as detected by a NanoDrop Spectrophotometer (Nano-100, Aosheng Instrument Co. Ltd., Hangzhou, China). For Illumina sequencing, we used PicoGreen ^® DNA Quantitation Assay (Thermo Fisher Scientific Inc., Waltham, MA, United States) and Qubit Fluorometric Quantitation (Thermo Fisher Scientific Inc., Waltham, MA, United States) to quantitate DNA amounts.

For bacterial and archaeal 16S rRNA genes, the V4 region was amplified with the primer pair 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACT ACHVGGGTWTCTAAT-3′) (Caporaso et al., 2011); for ITS genes, the ITS2 region was amplified with the primer pair gITS7F (5′-GTGARTCATCGARTCTTTG-3′) and ITS4R (5′-TCCTCCGCTTATTGATATGC-3′) (Ihrmark et al., 2012), combined with self-designed barcodes to distinguish samples. The details of PCR amplification, amplicon purification, library preparation, Illumina MiSeq sequencing, and sequence processing were described previously (Feng et al., 2017; Wei et al., 2018). Sequence processing was conducted on an in-house pipeline¹ integrated with the necessary bioinformatics tools (Feng et al., 2017; Wei et al., 2018). UPARSE (Edgar, 2013) was used to remove chimeras and classify the sequences into operational taxonomy units (OTUs) at a similarity of 97% with singletons being discarded. A randomly resampled OTU table was obtained to normalize total reads with a resample size of 24,084 reads for 16S rRNA gene and 10,253 reads for ITS gene based on the minimum sequence number per sample. The sequencing data are available in National Microbiology Data Center² with project number NMDC10017898.

The biodiversity indices for α-diversity and β-diversity were evaluated for both prokaryotic and fungal groups. The Shannon index, inverse Simpson index, observed richness, and Faith’s phylogenetic diversity were estimated to indicate the microbial species diversity. Based on weighted UniFrac distance matrix, principal coordinate analysis (PCoA) was used to show the difference among the prokaryotic or fungal communities from upper and lower soils under ambient and elevated temperature treatments, and the significance for warming and strata contributions to the microbial community structure were analyzed by permutational multivariate analysis of variance (PERMANOVA). Moreover, linear discriminant analysis effect size (LEfSe) was used to find the specialists for the different communities at the OTU level with the criterion of LDA score of 3.0 and significance level of 0.05. Further, the relationship between biodiversity and environmental variables were evaluated with Mantel test and constrained canonical analysis (CCA), and the explained proportion of variance by environmental groups was demonstrated with variance partition analysis (VPA). In addition, one-way ANOVA was used to test the differences of multiple groups. All biodiversity and statistical analyses were conducted with an online in-house pipeline (see text footnote 1; Feng et al., 2017).

In order to quantify the underlying mechanisms of microbial community assembly, we used a newly developed method, called Infer Community Assembly Mechanisms by Phylogenetic-bin-based null model (iCAMP) (Ning et al., 2020). iCAMP assigns the microbial OTUs into phylogenetic-closed groups (bins) and then quantifies the contribution of each ecological process to microbial community assembly, including homogeneous selection (HoS), heterogeneous selection (HeS), dispersal limitation (DL), homogenizing dispersal (HD), and undominated (DR). Afterward, the proportion of the ecological process between different samples can be compared based on the bootstrapped methods, and the contribution of each phylogenetic bin can be identified as well. In addition, another approach, normalized stochasticity ratio (NST), was used to determine the relative importance of determinism and stochasticity for microbial community assembly (Ning et al., 2019). The analyses were conducted with the “iCAMP” (v. 1.3.4) and “NST” (v. 3.0.6) packages in R (v. 4.1.1).

As our previous study showed (Wang et al., 2019), after 3 years of warming treatment, the soil moisture had increased. Warming also significantly increased mean plant height and plant biomass. Warming increased soil mean annual temperature 1.40 and 1.16°C at upper and lower soil strata, respectively; while it made no difference between upper and lower soil strata (Supplementary Figure 1). However, it did not alter other soil properties (Supplementary Figure 1) or alpha diversities of the prokaryotic (Supplementary Figure 2) and fungal communities (Supplementary Figure 3). Meanwhile, there were significant differences in soil properties between the upper and lower soil strata (Supplementary Figure 1), as well as alpha diversity of the prokaryotic (Supplementary Figure 2) and fungal communities (Supplementary Figure 3). Soil ammonium and nitrate content, total and available nitrogen content, and organic matter were significantly higher in the upper strata, while pH was the only property that was significantly lower in the upper strata. Alpha diversities of the prokaryotic community were higher in the upper soil, while those of the fungal community were higher in the lower soil strata.

The prokaryotic and fungal community compositions were significantly different between soil strata under unwarming treatment as revealed by PERMANOVA and PCoA (Table 1 and Supplementary Figure 4). However, microbial community compositions of the two soil strata were similar under warming treatment.

Based on the phylogenetically informed classification results from RDP classifier, linear discriminant analysis (LDA) effect size (LEfSe) was utilized to identify potential specialist taxa that were different between the two strata with warming. LEfSe identified 16 prokaryotic taxa above order level from unwarming lower soil, nine from unwarming upper soil, 10 from warming upper soil (Figure 1A and Supplementary Figure 5A). In addition, LEfSe identified 31 fungal taxa above family level from unwarming lower soil, one from unwarming upper soil, two from warming lower soil, and four from warming upper soil (Figure 1B and Supplementary Figure 5B). For the prokaryotic community, Gemmatimonadota and Actinobacteriota were more enriched in the unwarming lower soil, Alphaproteobacteria and Bacteroidota in the unwarming upper soil, while Firmicutes was more enriched in the warming upper soil (Figure 1A and Supplementary Figures 5A, 6A). Within the fungal community, Mortierellomycota and Eurotiomycetes were more enriched in the unwarming lower soil, Sordariomycetes in the unwarming upper soil, while Candida was more enriched in the warming lower soil and Ustilaginomycetes in the warming upper soil (Figure 1B and Supplementary Figures 5B, 6B). The above results indicated that more specialists were only observed under ambient control conditions than the warming treatment, suggesting that warming would reduce the differences between soil strata for prokaryotic and fungal specialists.

Mantel test results showed that the prokaryotic community was significantly related to all soil properties except ammonium and nitrate content, while the fungal community was only significantly related to total nitrogen content, but with a low r value (Table 2). The CCA model was significant for the prokaryotic community when the significant soil measurements were included as independent variables in explaining the variation of the prokaryotic community composition (p = 0.001; Supplementary Figures 7A,C), but was not significant for fungal community composition (p = 0.067; Supplementary Figures 7B,C). VPA results showed that soil properties and their interactions explained 28.16% of the variation in prokaryotic community composition (Supplementary Figure 7A), and 25.12% of the variation in fungal community composition (Supplementary Figure 7B). These results suggested that the prokaryotic community composition significantly correlated with most soil properties in this Tibetan alpine grassland, while fungal community composition did not. Therefore, the community structure of prokaryotes was more sensitive to environmental variables than fungi. However, soil properties explained only a small portion of the microbial community variation.

In order to explore the relative importance of determinism and stochasticity for microbial community assembly, the null-model-based NST method estimated the normalized stochasticity ratio and found that stochastic processes dominated the microbial assembly for prokaryotes and fungi across the treatments (Supplementary Figure 8). Essentially, the taxonomic stochastic ratio decreased for the prokaryotic community under warming from 79.7–82.9 to 62.5–65.0%, and there was a slightly higher stochastic ratio for the lower soil strata communities compared to upper soil samples (Supplementary Figure 8A). Fungal communities showed higher stochastic ratio with phylogenetic NST as compared with taxonomic NST, and warming also increased the deterministic ratio (<50%, Supplementary Figure 8B). These results suggested both prokaryotic and fungal communities showed less stochasticity under warming treatment and in the upper soils of the alpine grassland.

Furthermore, the iCAMP method assigned the deterministic and stochastic processes into five specific ecological processes (Figure 2). Homogeneous selection (58.1–61.9%) and dispersal limitation (30.1–42.5%) dominated the main processes for prokaryotic and fungal communities, respectively, for both warming treatment and strata profiles (Figures 2A,B), indicating higher deterministic components for prokaryotic rather than fungal community assembly. Moreover, warming significantly increased the proportion of the deterministic component (heterogeneous selection and homogenous selection) both for prokaryotic (58.6–59.1 vs. 61.8–62.8%) and fungal (26.4–29.6 vs. 33.2–34.4%) communities within the same strata (Student’s t test; p < 0.001), showing deterministic tendency under climate warming scenario (Figures 2C,D), especially for the homogenous selection process within the prokaryotic (58.1–58.6 vs. 61.3–61.9%) and fungal communities (26–29.5 vs. 31–32.7%). As for the stochastic component, the proportion decreased in fungi (70.4–73.6 vs. 65.6–66.8% vs.) under warming treatment. Considering strata profiles, the ecological processes showed a distinct tendency within the prokaryotic and fungal communities. For example, the proportion of dispersal limitation increased (31.1–31.9 vs. 36.2–42.5%) and the proportion of homogenous selection decreased (29.5–32.7 vs. 26–31%) within lower soils as compared to upper soils for the fungal community, while the prokaryotic community showed the opposite tendency. More importantly, the difference in stochasticity between upper and lower fungal communities decreased under the elevated temperature condition (65.6 vs. 66.8%) as compared to ambient temperature condition (70.4 vs. 73.6%), and dispersal limitation displayed a similar trend (Figure 2). Thus, the prokaryotic and fungal communities showed different assembly mechanisms between soil strata for ambient and elevated temperature conditions.

To further understand the underlying mechanisms of such distinct ecological process patterns, we looked into the ecological processes of the dominant groups that may make the greatest contributions to community variance, and the relationships between ecological processes and environmental parameters. First, the representative ecological processes were identified for each bin, which were grouped in closed phylogenetic relations, and the change of ecological process for each bin under different conditions may explain the differences in community assembly for prokaryotes and fungi at a finer scale (Figure 3 and Supplementary Figure 9). The three bins with higher relative abundance within the whole community for prokaryotes were bin27 (11.5%), bin115 (5.8%), and bin72 (4.8%), of which all dominant processes were assigned to homogeneous selection between the two strata and temperature conditions (Supplementary Figure 9). However, the representative processes of the top three fungal bins, i.e., bin12 (22.0%), bin4 (10.7%), and bin24 (6.0%), with high relative abundance changed between the different strata and temperature conditions (Figure 3). The ecological process for bin12 was undominated for unwarming lower soils, while all others were homogeneous selections. Similarly, warming changed the assembly mechanism of bin24, from undominated to dispersal limitation. These variations in status for different bins may be associated with assembly of the whole microbial community. Considering the influence of soil variables on the microbial community, Mantel test was used to correlate the soil heterogeneity with microbial ecological processes (Supplementary Tables 2, 3). Soil nutrient elements (e.g., total nitrogen, available nitrogen, and organic matter) were significantly correlated to homogeneous selection of prokaryotes under ambient temperature condition, and dispersal limitation under elevated temperature condition tended to be associated with nutrient elements as well (Supplementary Table 2). Different from the prokaryotic community, the ecological processes of fungi showed less correlation to environmental variables, but the dispersal limitation was slightly correlated to ammonia nitrogen and organic matter for lower soils at ambient (r_mantel = -0.341, p = 0.0419) and elevated (r_mantel = 0.421, p = 0.0105) temperature conditions, respectively (Supplementary Table 3). The differences in soil nutrient elements between upper and lower soils may be associated with microbial assembly processes, especially the difference for ambient and elevated temperature condition.