Deterministic processes dominate soil microbial community assembly in subalpine coniferous forests on the Loess Plateau

Site and sampling

A total of 23 soil plots were sampled (Table S1, Fig. S1) in August 2016 and August 2017. The sites were selected because their vegetation was subalpine mountain coniferous forests and they were located between 1,900 m and 3,055 m above mean sea level (amsl). The study area has a warm temperate continental monsoon climate, and mostly cinnamon soil.

This study focused on response patterns along environmental gradients rather than exploring differences among treatment groups. Thus, we sampled along three altitudinal gradients without replicates. Previous studies have shown that for continuous environmental drivers, gradient designs further allow for better extrapolation, characterization of (nonlinear) response functions, and, consequently, quantitative outputs better suited for ecological models than replicated designs (Cottingham, Lennon & Brown, 2005).

To avoid the interference of vegetation factors, we sampled plots in the single vegetation type (i.e., Larix principis-rupprechtii forests). These sites located on subalpine ecological environments possess pronounced climatic gradients and climosequences within short distances, with a high level of environmental heterogeneity (Siles & Margesin, 2017). Therefore, the sites with different altitudinal gradients corresponded to different environmental contexts and different environment characteristics.

Eight plots were sampled from the Wutai Mountain site (WT), which ranges from 1,900 m and 3,055 m amsl. Ten plots were sampled from the Pangquangou Natural Reserve site (PQG), ranging from 1,950 m to 2,650 m amsl. The last, five plots were sampled from the Luya Mountain site (LY), which ranges from 2,000 m and 2,400 m amsl. The details of each sample plots were added in the Supplemental Files (Table S1). At each sampling site, a 1 m × 1 m sampling plot was established in situ along the elevation gradient. Five soil cores at a depth of 15 cm were taken at each sampling plot, and then combined to form a single independent soil sample. Then, the soil samples were sealed in plastic bags and refrigerated, immediately transported to the laboratory and sieved using a 2 mm mesh. The soil samples were then stored at −80 °C until further analysis.

The soil samples were subsampled for molecular analysis and the DNA from of 1 g of soil was extracted using an E.Z.N.A.@ Soil DNA Kit (Omega Bio-tek, Inc., Norcross, GA, USA). The quality and quantity of the DNA extracts were measured using an Infinite 200 PRO plate reader (TECAN, Männedorf, Switzerland). The DNA purity was assessed based on the A260/A280 absorbance ratios, and only DNA extracts with absorbance ratios of 1.8∼2.0 were used for further analyses. Three DNA samples were extracted, from each soil sample, which were then combined and sequenced at Shanghai Personal Biotechnology Co., Ltd. on an Illumina MiSeq sequencing platform based on the bacterial v3–v4 hypervariable region using bacterial 16S universal primers (341F 5′-ACTCCTACGAGGAGCA- 3′ and 805R 5′-TTACCGCGGCTGCTGGCAC- 3′) (Tripathi et al., 2018).

Bioinformatics analysis

The sequencing data were analyzed using the QIIME pipeline (v1.8.0, http://qiime.org/) (Caporaso et al., 2010). The filtered sequence alignments were denoised by DeNoiser (Reeder & Knight, 2010) and then screened for chimeras using UCHIME (Edgar et al., 2011). The Archaea and unknown sequences were removed. The sequences were clustered into operational taxonomic units (OTUs) at a 97% similarity level using the average neighbor method and taxonomy was blast to SILVA database by k-mer searching using MOTHUR (Pruesse et al., 2007). The OTU table was rarefied to 4,020 sequences per sample. Ten independent maximum-likelihood phylogenetic trees based on Jukes–Cantor distance were then constructed using FastTree2 (Price, Dehal & Arkin, 2009) after the removal of gaps and hypervariable regions using a Lane mask supplied by QIIME to support phylogenetic diversity calculations.

Environmental variables

In the laboratory, soil total carbon (TC), total nitrogen (TN), and total sulfur (TS) were measured using an elemental analyzer (Vario EL/ MACRO cube, Elementar, Hanau, Germany); nitrate nitrogen (NO₃⁻_N), ammonium nitrogen (NH_4⁺_N), and nitrite nitrogen (NO₂⁻_N) were measured by an Automated Discrete Analysis Instrument (CleverChem 380; DeChem-Tech, Hamburg, Germany). After shaking the soil: water suspension (1:2.5 mass/volume) for 30 mins, the soil pH was measured using a pH meter (Hl 3221, Italy). The soil organic carbon in each soil sample was measured using the potassium dichromate volumetric method (Nelson & Sommers, 1982).

Null model analysis

A null model was constructed to account for changes in β-diversity while controlling for stochastic variation and associated changes in α-diversity (i.e., local species richness; 999 iterations) (Chase et al., 2011). We considered the null deviation as the relative difference between the observed β-diversity and the null-model β-diversity (Tucker et al., 2016). As such, null deviation values may represent communities that are more similar than expected by chance (a negative null deviation value), less similar than expected by chance (a positive null deviation value), or close to the chance expectation (values near zero) (Tucker et al., 2016).

Phylogenetic analysis

Our study used phylogenetic turnover between communities to infer ecological processes (Stegen et al., 2015). To quantify phylogenetic turnover between communities, we used the between community mean-nearest-taxon-distance (βMNTD) metric. βMNTD was calculated in R (R Core Team, 2018) ‘comdistnt’ (abundance.weighted = TRUE; package “picante”). Then, we evaluated β-Nearest Taxon Index (βNTI), which expresses the difference between observed βMNTD and the mean of the null distribution in units of standard deviations (Stegen et al., 2013).

In addition, to distinguish more details in the assembly processes, we used the Raup–Crick metric (Chase et al., 2011), extended to incorporate species’ relative abundances; referred to as RC_bray. The R script of RC_bray can be found at https://github.com/stegen/Stegen_etal_ISME_2013.

In a given community, we estimated the relative influence of variable selection or homogeneous selection as the fraction of their comparisons with βNTI > +2 or βNTI < −2, respectively. We regard the fraction of the between community comparisons with |βNTI| < 2 and RC_bray > +0.95 as dispersal limitation, while |βNTI| < 2 and RC_bray < −0.95 is considered homogenizing dispersal (Diniandreote et al., 2015; Stegen et al., 2013; Stegen et al., 2015).

Network analysis

The co-occurrence network was constructed based on the Spearman correlation matrix offered in the ‘psych’ package in R. In this network, the nodes represent OTUs and the edges that connect these nodes represent correlations between OTUs. Only those connections with correlation coefficients >0.6 and P < 0.05 were used in the network. Thus, positive correlations indicate co-occurring OTUs based on abundances, whereas negative correlations indicate that the OTUs are mutually exclusive (Barberán et al., 2012). P-values were false discovery rate (FDR) adjusted to control for the analysis (FDR < 0.05). The network analysis was completed using the ‘igraph’ package in R.

Statistical analysis

All statistical analyses were performed in the R environment using the ‘vegan’, ‘ggplot2’, ‘ggpubr’, and ‘corrplot’ packages. A Venn diagram was used to visualize the shared OTUs among the sites. A correlation matrix graph was used to demonstrate the correlation between soil physicochemical factors and was constructed using the ‘corrplot’ packages in R. Multivariate regression tree analysis (MRT) was used to explain the relationship between bacterial α-diversity estimates and environmental variables in a visualized tree, and diversity indices were normalized to the same mean before performing MRT analysis (Ge et al., 2008). Based on the longest gradient lengths from the results of detrended correspondence analysis (DCA), we selected redundancy analysis (RDA) to quantify the effects of environmental variables on microbial community composition (Mo et al., 2018). Forward selection of PCNM variables based on permutation tests was chosen to identify two of the 23 extracted PCNM variables that significantly (P < 0.05) explained the spatial structure. The PCNM eigenfunctions, which represent the ‘spectral decomposition of the spatial relationship across sampling locations’, can be considered as the spatial variables in the ordination-based analysis. The contributions of environmental filtering and the space variable (PCNM) to the variation in bacterial community composition were calculated by using variance partitioning analysis (VPA) (CANOCO for Windows Version 5.0). The mantel test was performed in the R environment using the ‘vegan’ packages.

Physicochemical properties of the soils from the different sites

The soil physicochemical properties varied across the different sampling sites (Fig. 1). Briefly, the contents of ammonium nitrogen and nitrite nitrogen were the highest at LY sites (36.91 and 0.16 mg kg⁻¹, respectively), and were lowest at WT sites (17.41 and 0.04 mg kg⁻¹, respectively). The contents of nitrate nitrogen (6.45 mg kg⁻¹), SOC (70.29 mg g⁻¹), TC (6.4%), and TN (0.51%) were the highest at WT sites, and were the lowest at LY sites.

TN was significantly positively correlated with TC and SOC (P < 0.05) and significantly negatively correlated with pH value (P < 0.05; Fig. 2). TC and pH showed a significant negative correlation (P < 0.05). SOC was significantly positively correlated with nitrate nitrogen (P < 0.0.5) and significantly negatively correlated with nitrite nitrogen (P < 0.05). This indicated that the sites sampled had different environment characteristics.

Dynamics of bacterial community composition and diversity

A total of 4,258 OTUs were identified from 1,062,241 high-quality sequences recovered from 23 soil samples. Good’s coverage ranged from 95.19% to 99.75%, indicating that the identified sequences were representative of most of the bacterial sequences in the collected soil samples. Rarefaction curve analyses, which generally yielded asymptotic curves, indicated that the numbers of sampling plots were enough. Detailed information of the sequencing results is provided in Table S2.

The soil microbial community composition varied across the different sampling sites (Fig. 3). There were 15 bacterial phyla with relative abundances of more than 0.01% (Fig. 3A). As shown in the Venn diagram, 869 bacterial shared OTUs were observed in all sampling sites. There were 46 bacterial phyla identified (Fig. 3B). The abundance of Proteobacteria at all sites was the highest (mean relative abundance = 30.59%), and followed by Acidobacteria (19.63%), Actinobacteria (16.51%) and Chloroflexi (13.22%). Briefly, the mean relative abundance of Proteobacteria was the most at PQG (34.39%), and that of Actinobacteria was the highest at LY (26.29%). The mean relative abundances of Acidobacteria (28.68%) and Chloroflexi (16.09%) were the highest at WT. There were 31 bacterial families with relative abundances of more than 0.01% (Fig. 3C). Based on the clustering graph, the sampling plots of each of the sites roughly clustered together (Fig. 3D). The community α-diversity indices varied at the different sites (Fig. 4). Briefly, the phylogenetic diversity (pd) and the number of observed species (sobs) were the highest at WT sites (P < 0.05). There was no significant difference in the ACE index, Chao index, Shannon index and Simpson index at the different sites (P > 0.05). This indicated that the sites sampled had different soil microbial community structure.

Effects of environmental factors on microbiome dynamics

Based on the results of the DCA (axis length = 1.02), we used RDA to identify the abiotic environmental drivers that influenced bacterial community composition (Fig. 5; permutation test, P < 0.01). The results demonstrated that Proteobacteria, Bacteroidetes, and Cyanobacteria were mainly driven by pH, while SOC, TC, and TN were the main abiotic drivers of Parcubacteria and Planctomycetes.

In the MRT analysis (Fig. 6), we observed that the diversity indices (normalized) were mainly split by SOC, explaining 36.75% in the first spilt. The correlation analysis showed similar results: SOC was significantly correlated with bacterial communities at the phylum level (e.g., Proteobacteria, Bacteroidetes, and Chloroflexi). Given its contribution to explaining community distribution patterns, SOC was further used as a descriptor for the environmental gradients.

The variation partitioning analysis showed that environmental variables (20.3%) explained more variation of microbial community structure than spatial variables (1.9%). This suggested that both deterministic and stochastic processes were involved in the assembly of microbial communities, and that deterministic processes were dominant. The unexplained variable was 78.6% (Fig. 7).

Nonrandom co-occurrence patterns of the microbial community

Network analysis was applied to explore the interspecific relationship patterns in the microbial communities (Fig. 8). Compared with the LY- and WT- network, the PQG-network exhibited more edges (87), more vertices (40), more modularity (0.691), higher average degree (4.35) and average clustering coefficients (0.858), but less the numbers of modules (6) (Table S3). Strong positive correlations were observed at all sites, while negative correlations were rare. The size of the nodes corresponds to betweenness centralization values.

The bacterial community assembly processes

According to the null model analysis, our results demonstrated that the null deviation values varied at different sites (ranging from 0.29 to 0.57; Fig. 9A). The bacterial communities at WT deviated significantly from the null expected value (relative null deviation = 0.45) and were greater than that at LY site and PQG site (relative null deviation = 0.32 and 0.34, respectively) (P < 0.05).

Most importantly, we observed that the microbial community was more greatly shaped by variable selection (βNTI > +2) (Fig. 9B). From LY to WT, we observed a gradual increase in the relative role of deterministic processes compared to stochastic processes (Fig. 9C). Based on the regression analysis of the environmental variables with assembly process parameters, we found that SOC had a great influence on community assembly processes (Fig. 9D). The mantel test between βNTI and SOC matrices indicated the similar conclusion (P < 0.05, R = 0.509).