Proximate grassland and shrub-encroached sites show dramatic restructuring of soil bacterial communities

Site description and sample collection

The study area was selected in a high-density shrub-encroached grassland (42°57′N, 112°43′E; 1,208 m; Fig. S1), located in Inner Mongolia, China. The average annual temperature is 5.1 °C and the mean precipitation is 195 mm in this region (Chen et al., 2015). The dominated grass is Cleistogenes songorica across the region, but Caragana microphylla is encroaching (Chen et al., 2015). Soil samples were collected on the 10th of August, 2016. We identified ten shrub-encroached sample plots to include in this study. The selected sites were more than 500 m away from each other. At each site (10 × 10 m), the encroachment soils were sampled under five shrub patches (the nearest to the four vertices and the center of a plot) with 0–10 cm depth and mixed as one sample. The control non-encroached soils were collected one m away from the five shrub canopies with 0–10 cm depth and mixed as one sample (Fig. S1). In total, 10 from control grassland soils and 10 from adjacent encroached soils were collected for further study. The soils were fully mixed and sieved, and then transported refrigerated to the lab within 24 h. The soils were divided into two parts: one part was stored at 4 °C for biogeochemical analysis and the other was stored at −20 °C for DNA extraction.

Sample pretreatment

Measurement of soil properties, DNA extraction, and amplicon library preparation are described in the Supplemental Information.

Processing of sequence data

The raw data were processed by QIIME (v.1.9.0; Caporaso et al., 2010). The sequences were clustered into operational taxonomic units (OTUs; 97% identity) with UCLUST (Edgar, 2010). Chimeric and singleton OTUs were removed prior to downstream analysis. The default setting was used to select the representative sequence (i.e., most abundant sequence) for each OTU, which was assigned taxonomic annotations using the UCLUST (Edgar, 2010) and aligned by PyNAST (Caporaso et al., 2010). To normalize for sampling depth, random subsets of 26,000 reads per sample (the lowest sequence read depth across the study) were used to calculate bacterial alpha- and beta-diversities.

Statistical analysis

Phylogenetic diversity (PD) was estimated by Faith’s index (Faith, 1992). Pairwise t-test was performed to show differences in relative abundance of dominant bacterial phyla and alpha-diversity. Pearson correlation was used to test relationships between bacterial alpha-diversity and soil properties. Linear discriminant analysis effect size (LEfSe) was used to identify bacterial taxa that differed significantly between treatments (default setting; Segata et al., 2011). Non-metric multidimensional scaling and Analysis of Similarity (ANOSIM; permutations = 999) were performed to distinguish the differences in bacterial community composition between treatments by using the vegan package (v.2.0-2) in R software. The correlation between variables (i.e., soil properties and spatial distance) and soil bacterial community composition were analyzed by Mantel tests (permutations = 999). Multicollinearity of soil properties was tested by the variance inflation factor (VIF; Zuur, Leno & Elphick, 2010), and those properties with the VIF values < 3 were selected for canonical correspondence analysis (CCA).

The nearest taxon index (NTI) and beta nearest taxon index (betaNTI) were performed using the picante package (Purcell et al., 2007) and Phylocom 4.2 (Hardy, 2008), respectively, to analyze soil bacterial phylogenetic structure. The NTI measures the mean nearest taxon distance among individuals to estimate the phylogenetic dispersion of the community (Webb, 2000). More positive or negative NTI values indicate phylogenetic clustering or overdispersion, respectively (Webb, 2000). BetaNTI values between −2 and 2 suggested stochastic process (neutral assembly) while the values above 2 or below −2 indicated deterministic processes (niche assembly, Stegen et al., 2012). Co-occurrence networks were generated in R using the “WGCNA” package (Langfelder & Horvath, 2012). We adjusted all P-values (cutoff as 0.001) by using the Benjamini and Hochberg false discovery rate for multiple testing (Benjamini, Krieger & Yekutieli, 2006). The network nodes defined as network hubs (z-score > 2.5; c-score > 0.6), module hubs (z-score > 2.5; c-score < 0.6), connectors (z-score < 2.5; c-score > 0.6), and peripherals (z-score < 2.5; c-score < 0.6) referring to their roles in network structure (Poudel et al., 2016). Network hubs are those OTUs that are highly connected both in general and within a module. Module hubs and connectors are OTUs that are highly connected only within a module and only link modules, respectively. Peripherals are defined as OTUs that have few links to other species. The bacterial metabolic function was predicted by phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) according to KEGG database (Langille et al., 2013).

Soil chemistry

Compared to non-encroached grassland soils, shrub-encroached soils were associated with higher content of NO₃⁻, total nitrogen (TN), total carbon (TC), dissolved organic carbon (DOC), soil organic carbon (SOC), and total phosphorus (Table S1). However, shrub encroachment showed little effect on other soil properties, such as soil pH, NH₄⁺ content and SM relative to control in this study (Table S1).

Bacterial alpha-diversity

A total of 966,631 quality bacterial sequences was obtained with 26,037–68,261 (mean 48,332) sequences per sample. In this study, bacterial alpha-diversity included OTU richness, Shannon index, evenness, and PD, which was calculated by randomly selected subsets of 26,000 reads per sample. Generally, encroached sites had significant higher alpha-diversity relative to grassland sites (Fig. 1). Bacterial OTU richness was positively correlated with NO₃⁻, DOC, TC, TP, and SOC; PD was positively correlated with NO₃⁻, DOC, TC, and SOC; the Shannon index was positively correlated with NO₃⁻, TC, TP, and SOC; evenness was positively correlated with NO₃⁻, TN, TC, and SOC (Table 1).

Bacterial community structure

The dominant soil bacterial phyla (i.e., relative abundance > 1%) across all samples were Actinobacteria (27.3%), Acidobacteria (23.1%), Proteobacteria (23.0%), Chloroflexi (6.0%), Planctomycetes (4.7%), Gemmatimonadetes (2.8%), Firmicutes (2.7%), Bacteroidetes (2.6%), and Nitrospirae (2.4%) (Fig. S2). Compared to control grassland soils, the relative abundance of Proteobacteria showed significantly lower in encroached sites (Fig. S3). Compared to control, encroachment was associated with higher relative abundance of Chloroflexi and Nitrospirae (Fig. S3). LEfSe analysis showed that bacteria in one phylum (i.e., Proteobacteria), five classes (i.e., Acidobacteriia, ML635J_21, vadinHA49, Solibacteres, and Gammaproteobacteria) and 15 orders (i.e., Acidobacteriales, Solibacterales, Planctomycetales, Caulobacterales, Rhodospirillales, Burkholderiales, etc) were significantly more abundant in control soils. Bacteria from two phyla (i.e., Nitrospirae and Armatimonadetes), two classes (i.e., Nitrospira and Chloroflexi) and 11 orders (i.e., Gaiellales, Roseiflexales, Nitrospirales, Syntrophobacterales, Desulfovibrionales, etc) were significantly more abundant in shrub-encroached soils (Fig. 2).

Significant differences in soil bacterial community compositions were found between shrub-encroached and grassland sites (ANOSIM: P = 0.001; Fig. 3). The NTI values showed positive (i.e., >0; P = 0.001) for all samples, indicating that bacterial phylogenetic structure showed clustering in both encroached and control soils (Fig. 4). Almost all betaNTI scores for bacterial communities were below −2, which suggested that deterministic assembly dominated soil bacterial community dynamics in both grassland and shrub-encroached soils (Fig. 4). A correlation network was built at bacterial genus level. There was a larger proportion of positive than negative correlations between genera in soils (Fig. S4A). Compared to grassland soils, shrub-encroached soil showed higher proportion of correlation network hubs (Fig. S4B), suggesting that bacterial community in shrub-encroached soils might be more interconnected than grassland soils.

Mantel tests demonstrated that soil bacterial community composition showed significant correlation with soil pH, SM, NO₃⁻, DOC, TC, TP, and SOC (Table 2; P < 0.05 in all cases). Among these variables, SOC content (P = 0.002) had the strongest correlation with soil bacterial community composition. However, spatial distance showed little correlation with bacterial community composition (P = 0.181; Table 2). CCA further demonstrated that SOC was the primary driver affecting soil bacterial community composition (Fig. S5).

The predicted metabolic function

The metabolic function of bacterial community was predicted by PICRUSt. A total of 328 predicted functional genes were detected in this study. More than 89% of total sequences belonged to categories of metabolism (52.2%), genetic information processing (15.8%), environmental information processing (13.3%), and organismal systems (8.35%) in soils, according to the KEGG database. Compared to controls, shrub encroachment was associated with significant differences in potential functions of the soil bacterial community (Fig. S6). Metabolism of cofactors and vitamins, energy metabolism, glycan biosynthesis and metabolism, enzyme families, and nucleotide metabolism were enriched in grassland soil, while xenobiotics biodegradation and metabolism, lipid metabolism, metabolism of terpenoids and polyketides, amino acid metabolism, and carbohydrate metabolism were enriched in shrub-encroached soils (Fig. 5). The relative abundances of sequences associated with cell motility, environmental adaptation, signal transduction, and protein folding, sorting and degradation were enriched in grassland soils (Fig. 5). The sequences related to cell growth and death, transport and catabolism, nervous system, membrane transport, and transcription were enriched in shrub-encroached soils (Fig. 5).