The effects of afforestation on soil bacterial communities in temperate grassland are modulated by soil chemical properties

Study sites and sampling

The study site was located 6 km west of the town of Honghua’erji in an artificial forest produced by seedling afforestation in a large area of thin grassland (savanna, N48.257443° E119.996448°). The annual precipitation is approximately 390 mm, and the mean annual temperature is roughly −2.4 °C. The elevation is approximately 740–1,100 m. The forest coverage is as high as 69.8% in the forest regions of the nature reserve, while the Honghua’erji grassland is a bare grassland ecosystem that is mainly composed of weedy species such as Stipa baicalensis, Festuca ovina, and Carex pediformis (Wen et al., 2002). The field experiments were approved by the Honghua’erji Nature Reserve (permit number 200/66150221). The tree ages of the Mongolian pines in this afforestation forest range from 27 to 33 years, with a diameter at breast height (DBH) of 21.21 ± 4.19 cm (12.09∼32.10 cm) and tree height of 12.86 ± 0.86 (11.2∼14.6 m).

At the study site, we collected soil pit samples from 20 quadrants, including 10 inside the forest (i.e., forest soils) and 10 in the adjacent grassland (i.e., grassland soils). Each surface area of quadrant was roughly 100 m². At each location, 1-kg soil samples were collected at a depth of 50 cm in permeable bags of nylon mesh from 10 quadrants located 10 to 100 m from each other (Fig. S1). This depth was chosen because the major microbiome composition and the major transition of soil physicochemical characteristics are located at a soil depth of approximately 50 cm (Fierer, Schimel & Holden, 2003). All equipment was sterilized by an autoclave. The soil samples were separated into two parts: one part was dried for quantification of soil physicochemical properties, and the other was used to quantify the microbiome. The latter samples were stored in RNAguardian stabilization solution (MBGEN Biosciences, Taipei, Taiwan) on ice immediately after collection until transfer to the laboratory, where the samples were stored at −20 °C before metagenomic DNA extraction.

Quantifying soil chemical properties

Organic carbon (C) was measured by the external-heat potassium dichromate oxidation method, and total nitrogen content (N) was measured by the Kjeldahl distillation method. The inductively coupled plasma (ICP) method was used with soils digested in a mixture of HF–HClO₄–HNO₃ to quantify the contents of the soil elements K⁺, P, Ca²⁺, Fe ³⁺, Mg²⁺, and Na⁺. Soil pH was determined using a Sartorius pH meter PB-10 (Germany) with a soil:water ratio of 1:2.5. All soil properties were determined and quantified by the State Key Laboratory of Vegetation and Environmental Change, Institute of Botany, Chinese Academy Sciences.

16S rRNA metagenome sequencing

Bacterial metagenomic DNA was extracted with an EZNA^® Soil DNA Kit (Qiagen, Valencia, CA, USA), and the concentration was adjusted to 50 ng ml⁻¹ via dilution. The metagenomic DNA was quantified using Qubit^® 2.0 (Invitrogen, Life Technologies, CA, USA). The primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′) were used to amplify the V3-V4 hypervariable 16S rRNA region (Mizrahi-Man, Davenport & Gilad, 2013), and the PCR products were used to construct a DNA library with the Roche GS FLX Titanium emPCR kit (Roche Applied Science). The DNA libraries were then sequenced by Sangon Biotech Co. (Shanghai, China) on an Illumina MiSeq 2X300. The sequencing procedures followed the manufacturer’s instructions.

Before analysis, the raw HTS data were cleaned by removing sequence fragments shorter than 200 bp or with missing barcodes or polyN or polyA/T in the Ribosomal Database Project (RDP) (Cole et al., 2007). We also discarded reads with PHRED quality scores <Q25 (Ewing & Green, 1998; Ewing et al., 1998). The Mothur package was used to remove non-prokaryotic sequences and de-noise and trim the sequences (Mothur ver. 1.30.1, http://www.mothur.org/wiki/MiSeq_SOP (Schloss et al., 2009). Chimeric sequences were removed using Uchime (Edgar et al., 2011). After quality filtering, we clustered sequences using a criterion of >97% sequence similarity as operational taxonomic units (OTUs) defined as representing the same species by the UPARSE pipeline (http://drive5.com/uparse/) (Edgar, 2013) and closed-reference OTU picking. We estimated the relative abundances (RAs) of the soil microbiome according to 16S rRNA metagenome sequencing. The rarefied OTU table was generated using Qiime (Caporaso et al., 2010) and deposited in Mendeley (doi: 10.17632/gjskh8wswz.1). All steps were conducted by Sangon Biotech Co. (Shanghai, China) according to their pipeline (https://www.sangon.com/services_ngs_metseq.html). Each OTU was annotated and classified according to the RDP classifier and SILVA database. The raw sequence data were deposited in NCBI GenBank under Bioproject PRJNA317430 (accession number: SAMN04607375).

Predictive functional profiling of bacterial communities

PICRUSt v. 1.0.0, a functional prediction tool for estimating shared gene content according to the corresponding bacterial phylogeny, was used to predict the functional potential of each sample (Langille et al., 2013). PICRUSt generates the composition of gene families for each metagenome using an extended ancestral-state reconstruction algorithm. The online version of PICRUSt implemented in Galaxy (http://huttenhower.sph.harvard.edu/galaxy/) was used to assist the algorithms. The quality-filtered sequences were assigned to an OTU table against the Greengenes v. 13.5 OTU database (DeSantis et al., 2006) for PICRUSt prediction implemented in QIIME v. 1.8.0 (Caporaso et al., 2010). Each OTU was normalized by its copy number. The functional contribution of each OTU member was reconstructed and predicted by mapping the 16S sequences to their nearest reference genome. ‘Virtual’ metagenomes with gene content abundance were then generated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog and Clusters of Orthologous Groups (COGs) databases.

Statistical analyses

All statistical analyses were conducted in R (R Core Team, 2013). This study included two types of measurements, the soil chemical properties and soil bacterial communities (i.e., relative abundances of bacterial OTUs), and one treatment, grassland afforestation (i.e., vegetation type: grassland vs forest). Due to inequality of variance (certain factors such as Na⁺ and Ca²⁺ deviated from equal variance in Levene’s test), we first assessed the differences in soil properties between vegetation types by the nonparametric Mann–Whitney U (MW) test using R. In addition, to determine whether the diversity of prokaryotic microbiome differs after afforestation, we calculated diversity indexes such as the Shannon-Wiener H, reciprocal Simpson’s index 1/D, species richness, and Pielou’s evenness for each of the soil bacterial communities and compared them between afforested and grassland sites by the Mann–Whitney U test. We also adopted the R package DESeq2 (Love, Huber & Anders, 2014) to identify OTUs that were differentially abundant between the afforested and grassland sites. A logistic regression model was further used to determine if each of the soil properties predicted the grassland-afforestation treatment. The likelihood ratio test (LRT) was performed to identify the best-fitting model of the logistic regression. To test the hypotheses of a stochastic process of bacterial assembly (random distribution model) or resource-governed assembly (niche-based mechanism), we used Rank-Abundance Dominance (RAD) analysis to display logarithmic relative species abundances against species rank order (McGill et al., 2007). Bray–Curtis distance-based redundancy analysis (dbRDA) implemented in the R package vegan (Dixon & Palmer, 2003) was used to examine the explanatory proportions of the treatment (vegetation type) for these two measurements (soil properties and microbiomes). Because soil chemical properties can also affect the soil bacterial communities (Stutter & Richards, 2012), we further performed a partial dbRDA to assess the effects of grassland afforestation on the soil bacterial communities conditioned by soil chemical properties and the effects of the soil chemical properties on the soil bacterial communities conditioned by vegetation type using the Bray-Curtis distance. The best model were chosen via forward selection by the ordistep function implemented in the R package vegan (Dixon & Palmer, 2003). Type II ANOVA was used to evaluate the fit of the model of each constraint factor. We further identified divergent bacteria (i.e., forest and grassland specialists) using the supermajority rule (>2/3 difference in abundance) with the assistance of the multinomial species classification method (CLAM) test (Chazdon et al., 2011). The soil elements that significantly explained the soil bacterial communities by partial dbRDA were then used as independent factors to predict the abundance of bacterial specialists by the generalized linear model (GLM). We also performed hierarchical cluster analysis (H-cluster) and discriminant analysis of principal components (DAPC) to illustrate the grouping patterns between afforested and grassland sites using the R package adegenet (Jombart, 2008). Because the ecophysiological functions (COGs and KEGG modules) of the soil microbiomes were also predicted, all tests were repeated by replacing the soil bacterial OTUs with the relative abundances of the COGs and KEGG modules. All R codes and input files have been deposited as Mendeley data (https://data.mendeley.com/datasets/gjskh8wswz/1).

Effects of afforestation on soil chemical properties

To characterize the soil chemical properties, the contents of eight elements (organic C, total N, P, K⁺, Ca²⁺, Mg²⁺, Fe³⁺, and Na⁺) and the pH value of each soil sample were determined (Table 1). Among these soil variables, the contents of C, P, Ca²⁺, Mg²⁺, and Fe³⁺ were significantly higher in the forest soils than in the grassland soils (P < 0.01), and the pH of the grassland soils was significantly higher than the pH of the forest soils (P = 0.031) (Fig. 1). To estimate the explanatory proportions of the afforestation effect on the variance of soil properties, we performed dbRDA using vegetation type (i.e., forest and grassland) as the independent factor. Vegetation type significantly explained 50.18% of the variance of soil properties (P = 0.002, Table 2 and Fig. S2), indicating that afforestation has changed the soil properties.

If the soil elements are altered after afforestation, the contents of soil elements will predict the vegetation type (i.e., afforestation effect). To test this hypothesis, we compared the simple logistic regression (SLR) model using a single soil element as the independent factor and the multivariate logistic regression (MLR) model using all soil elements as independent factors (i.e., the full model M1) to the empty (null) model M0. Most of the SLR models and M1 rejected the null model M0 (P < 0.005), except the SLR models with K (LRT: P = 0.084) or Na (LRT: P = 0.108) as the independent factor. We further compared the SLR models with each single soil element (C, N, P, Ca²⁺, Mg²⁺, Fe³⁺, and pH) as the independent factor with the MLR model using all C, N, P, Ca²⁺, Mg²⁺, Fe³⁺, and pH as independent factors (full model M2) by LRT to determine which individual factors might be related to afforestation. The SLR models with Mg or Fe as the independent factor could not be rejected by M2 (P = 0.278 and 0.130, respectively, Fig. 2). Although the other SLR models with C, N, P, Ca²⁺, or pH as the independent factor were rejected by M2, none of the independent factors in M2 could significantly predict the presence of forestation (Z < 10⁻⁶, P > 0.9999 in each term). These results indicated that the contents of Mg²⁺ (0.456 ± 0.099 and 0.313 ± 0.035 g. kg⁻¹ in forest and grassland soils, respectively) and Fe³⁺ (8.355 ±1.047 and 5.607 ± 1.308 g. kg⁻¹ in forest and grassland soils, respectively) could singly reflect the changes in soil properties due to afforestation (Table 1 and Fig. 2).

Afforestation effects on the relative abundances of soil bacteria

Overall, a total of 694,993 sequences were obtained after rarefaction analysis (314,923 for grassland and 380,070 for forest). Divergence of the soil microbiome was inferred by significant or marginal differences in the diversity indices (Shannon–Wiener index H, W = 26, P = 0.0753; Reciprocal Simpson’s index 1/D: W = 23, P = 0.0432; Pielou’s evenness J: W = 8.5, P = 0.0019, Table 3). However, no difference was estimated in the species richness of the soil microbiome between forest and grassland (W = 50, P = 1, Table 3). Hence, we used RA as an indicator to compare the microbiomes of the soil samples (Fig. S3). H-cluster and DAPC (five first PCs of PCA used, which conserved 70% of the variance of the bacterial RA) presented similar patterns of clear divergence of soil microbiomes between the forest and grassland (Fig. 3). A total of 350 OTUs with an adjusted P-value <0.01 exhibited twofold changes between forest and grassland vegetation (178 OTUs were higher in grassland, while 172 were higher in forest). The differentially abundant OTUs were only from the most abundant phyla (Proteobacteria, Acidobacteria, and Verrucomicrobia) and Planctomycetes. These results confirm that the soil bacterial composition has changed due to grassland afforestation despite no change in bacterial richness.

Because the close-reference OTU picking strategy discards samples that do not bin with 97% similarity, assignments of different OTUs were rarely found among different samples. We performed the CLAM test (Chazdon et al., 2011) to classify these OTUs as four types of bacteria using the supermajority (2/3) rule: too-rare bacteria, generalist bacteria, and grassland- and forest-specialist bacteria. In this classification, a large proportion of OTUs (95.46%) belonged to the too-rare bacteria, which accounted for 38.42% and 36.73% of the abundance of grassland and forest soil bacteria, respectively; only 2.98% of the OTUs were generalists, but they accounted for 35.55% and 35.88% of the abundance of grassland and forest soil bacteria, respectively. Only 1.56% of the OTUs belonged to specialists, of which 1061 (0.7%) and 1224 (0.8%) bacterial OTUs were identified as forest and grassland specialists, respectively. The RA of grassland-specialist bacteria was 22.84% in grassland soils and 2.84% in forest soils, while the RA of forest-specialist bacteria was 3.19% in grassland soils and 24.55% in forest soils.

We further used the RAs of the bacterial OTUs as the dependent response to access the impact of afforestation on the soil bacterial communities. When vegetation type was used as the categorical independent factor, 16.3% of the variance of bacterial RA was explained significantly (P = 0.001). However, the significant explanatory effect of vegetation type was lost after removing (conditioning on) soil properties (P = 0.466, Table 2). When soil properties were used as a constraint factor, 64.4% of the variance of bacterial RA was explained by soil properties, in which organic C (11.4% explanation), P (10.0%), Ca²⁺ (9.0%), Mg²⁺ (7.1%), and Fe³⁺ (10.5%) significantly fit the model by type II ANOVA (Table 2). Similar significant explanations by C (6.7%), P (9.5%), Ca²⁺ (8.8%), Mg²⁺ (6.7%), and Fe³⁺ (5.4%) were obtained when the effect of vegetation type was removed (i.e., partial dbRDA, Table 2). In addition, when conducting forward model selection, the bacterial RA was best explained by vegetation type, Ca²⁺, Mg²⁺, Na⁺ and pH (constrained proportion: 48.07%). However, although pH is important in shaping the soil microbiome at continental scales (Lauber et al., 2009), the pH did not vary among vegetation types (forest: 6.3–6.36, grassland: 6.13–6.35) in our fine-scale study. There was no collinearity of pH with other factors [variance inflation factor (VIF): 1.125]. Therefore, we do not discuss the importance of pH in our study. These analyses suggested that afforestation has changed the soil chemical properties, with changes in at least Ca²⁺ and Mg²⁺ in all analyses, thus indirectly affecting the soil bacterial communities (Figs. S2 and S3).

Effects of soil properties on the divergence of bacterial phyla between forest and grassland soils

To identify the divergent soil bacteria, the habitat-specialist bacterial phyla were classified under the supermajority rule (i.e., 2/3 majority). At the phylum level, five phyla were habitat specialists, including one forest specialist (Thaumarchaeota, phylum RA: 2.3%) and four grassland specialists (Chloroflexi, Fibrobacteres, Nitrospirae, and Parcubacteria, 9.3%). The forest-specialist phylum Thaumarchaeota belongs to Archaea, while the four grassland specialists are Eubacteria.

We further tested the effects of soil properties on the abundance of these habitat specialists by GLM. Because organic C, P, Ca²⁺, Mg²⁺, and Fe³⁺ significantly explained the RA of soil bacteria in dbRDA, these five soil variables were used as independent predictors in GLM. The variances of the bacterial abundance of the sampled soils were significantly or marginally significantly greater than the means, suggesting overdispersion of these responses (P = 0.056, 0.001, 0.051, 0.002, and 0.002; alpha = 52.675, 25.871, 1.065, 23.101, and 76.612 in Thaumarchaeota, Chloroflexi, Fibrobacteres, Nitrospirae, and Parcubacteria, respectively). Hence, we used the quasi-Poisson model for GLM, which led to the same coefficient estimates as the standard Poisson model but with adjustment of the dispersion parameter for overdispersion. Under the quasi-Poisson model, Fe³⁺ content was marginally correlated with the abundance of the forest specialist Thaumarchaeota; for the grassland specialists, the P content was correlated with the abundances of all four phyla (Table 4). In addition, C, Ca²⁺, and Mg²⁺ were significantly correlated with Chloroflexi, Fe³⁺ was significantly correlated with Nitrospirae, and Ca²⁺ and Mg²⁺ were marginally correlated with Parcubacteria (Table 4).

Afforestation effects on the soil ecophysiological functions predicted by soil microbiome

To understand the ecophysiological functions of the soil bacterial communities, we predicted their functional composition using 16S rRNA gene and databases of reference genomes. A total of 4,659 clusters of orthologous groups (COGs) and 306 KEGG modules (Level-3 KEGG orthology) were identified. Similar to the analyses for testing the effects of afforestation and soil properties on the soil bacterial communities, we used vegetation type and soil elements as predictors to test the explanatory proportion and significance of each predictor on the RAs of the COGs and KEGG modules. dbRDA indicated that 22.0% and 21.2% of the variation of COGs and KEGG modules was significantly explained by vegetation type, respectively, whereas the explanatory proportion decreased to 3.4% and 4.5% when conditioning the soil-property effect (Table 5). Soil properties explained 83.0% and 82.4% of the variation of the COGs and KEGG modules in dbRDA, respectively, in which C, P, Ca²⁺, Mg²⁺, Fe³⁺, and N significantly or marginally fit the model according to type II ANOVA (Table 5). When conditioned on vegetation type, the explanatory proportion decreased slightly to 64.4% and 65.8% for the COGs and KEGG modules, respectively, and the remaining significant fitting factors were P, Ca²⁺, and Fe³⁺ (Table 5).

Among all COGs, 44 and 21 were identified as forest- and grassland-predominant COG categories by CLAM analysis, respectively. Most of the forest-predominant categories are involved in post-translational regulation, such as ubiquitin or E3 ligase, while many categories in grassland are membrane proteins or are involved in cell transport. By contrast, among all KEGG modules, only two forest-predominant categories and no grassland-predominant categories were identified. Further testing of the correlation of soil elements with these specialist COGs and KEGG modules revealed that 50 of the 65 COGs and both KEGG modules were significantly correlated with at least one soil element under Poisson or quasi-Poisson regression. These significant correlations indicate that soil properties, particularly P, Ca²⁺, and Fe³⁺, account for the changes in ecophysiological functional due to grassland afforestation.