Variations in Soil Bacterial Community Diversity and Structures Among Different Revegetation Types in the Baishilazi Nature Reserve

The field study was conducted at the Baishilazi Nature Reserve, in the mountainous region of the eastern Liaoning Province, China (40°50′00′′-40°57′12′′N, 124°44′07′′-124°57′30′′E). The total area of the Baishilazi Nature Reserve encompasses 7,407 hm² and belongs to the mountain range of the Changbai Mountain. This area is characterized by a continental monsoon climate, with warm wet summers, long cold winters, and strong diurnal temperature variation. The mean annual temperature is 6.4°C, with a mean annual precipitation of 1,158 mm. The characteristics of the eight stands are listed in Table 1; all forest stands were older than 40 years.

In July, 2017, we collected soil samples from Juglans mandshurica (JM), Quercus mongolica (QM), mixed broadleaf forest (MB), Betula ermanii (BE), conifer-broadleaf forest (CB), Larix gmelinii (LG), Pinus koraiensis (PK), and Pinus tabuliformis (PT) stands after removal of the litter layer. A total of 24 soil samples (three plots of 20 m × 20 m as three independent replicates) were taken from the A horizon of each stand with the use of a soil auger (8 cm in diameter, 10 cm deep). To ensure the representativeness of soil samples in each stands, the strip sampling method was used. The samples of 15–20 sampling plots within one stand were mixed together to obtain one composite sample per stand. Identically, three real replicate samples stand and placed in cooled boxes. In the laboratory, the samples were sieved through a 2-mm mesh to remove roots and other debris and subsequently divided into two parts. One part was stored at -80°C for DNA extraction, while the other part was air-dried for soil physicochemical analyses. Figure 1 shows the characteristics of the soils from all eight stands.

Soil DNA was extracted from 0.5 g soil, using the Fast DNA SPIN extraction kits (MP Biomedicals, Santa Ana, CA, United States), according to the manufacturer’s instructions. The NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) was used to measure the quantity and quality of the extracted DNAs.

A quantitative PCR amplification of the bacterial 16S rRNA genes V3–V4 region was performed, using the forward primer 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and the reverse primer 806R (5′-GGACTACHVGGGTWTCTAAT-3′). Sample-specific 7-bp barcodes were incorporated into primers for multiplex sequencing. The PCR amplifications were performed in two steps. During the first step, each of three independent 25 μl mixtures per DNA sample contained 5 μl of Q5 reaction buffer (5×), 5 μl of Q5 High-Fidelity GC buffer (5×), 1 μl (10 μM) of each forward and reverse primer, 2 μl (2.5 mM) of dNTPs, 0.25 μl of Q5 High-Fidelity DNA polymerase (5 U/μl), 2 μl (40–50 ng) of DNA Template, and 8.75 μl of ddH2O. Cycling conditions were 98°C for 5 min (one cycle), 98°C for 15 s, 55°C for 30 s, 72°C for 30 s (25 cycles), followed by 72°C for 5 min. The PCR amplicons were purified with Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN) and quantified using the PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, United States). After the individual quantification step, amplicons were pooled at equal amounts, and pair-end 2 × 300 bp sequencing was performed using the Illumina MiSeq platform with the MiSeq Reagent Kit v3. We have uploaded all raw sequences to the NCBI Sequence Read Archive under submission number SUB4489099 and BioProject number PRJNA489351↗.

Sequence data analyses were mainly performed using the software packages QIIME and R (v3.2.0). Operational taxonomic units (OTU)-level alpha diversity indices, such as the Chao1 richness estimator, the abundance-based Coverage Estimator (ACE) metric, the Shannon index, and the Simpson index, were calculated using the OTU table in QIIME. The shared and unique OTUs among samples were used to generate Venn diagrams with the R software package. The linear discriminant analysis (LDA) effect size (LEfSe) algorithm method was used to detect the potential biomarkers based on a normalized relative abundance matrix. The heatmap representation of the relative abundance of bacterial OTUs among samples was built using R. Principal components analysis (PCA) was also conducted based on the genus-level compositional profiles.

Analysis of variance (ANOVA) was performed using SPSS 19.0 software. Soil physicochemical characteristics as well as bacterial total diversity and abundance were compared using LSD tests. Pearson’s correlation analysis was accustomed to estimate the correlations between soil characteristics and bacterial diversity indices. Canonical correspondence analysis (CCA), performed via Canoco 4.5, was used to evaluate the linkages between dominant bacterial groups related to soil environmental factors.

We found significant differences among stands in terms of soil total Carbon (TC) and total Nitrogen (TN) (Figure 1). Interestingly, both TC and TN were highest in JM, with 100.53 and 7.80 g/kg, respectively; while in PT, we measured only 31.76 g/kg and 2.48 g/kg, respectively. Available N (AN) followed the order of JM > MB > QM > BE > CB > LG > PK > PT. In all sites, the soil C/N ratio was below 25:1, with CB showing the highest value. Soil pH value ranged from 4.89 to 5.70; under QM, pH was the lowest with 4.89, followed by CB, while JM had the highest soil pH value. Significant differences were found for available P (AP) and total P (TP), with the highest values under PK and CB (Figure 1).

We obtained a total of 1,448,252 bacterial sequences, which remained after removing low quality sequences and chimeras (on average, 60,343 per sample). In total, 57,139 to 65,460 sequences per sample were obtained and classified into the domain level (Table 2); average sequence length was 330 bp. At the phylum, the remaining sequences were identified for 65,479 OTUs. A maximum of 2,845 OTUs were detected in CB; however, only 2,554 OTUs were obtained in the JM.

Venn diagrams were used to compare the bacterial communities under JM, QM, LG, CB, and PT, based on shared and unique OTUs among the samples. The number of shared OTUs among JM, QM, LG, CB, and PT was 1,187. A total of 392 unique OTUs were found in CB, while 344 unique OTUs were found in LG, 312 in PT, 199 in QM, and 407 in JM. In terms of shared OTUs, 290 were observed between LG and PT, and 502 between QM and JM. Under QM, the lowest number of unique OTUs was found (Supplementary Figure), while we detected most unique OTUs under JM. S1

The different soils showed high species richness, albeit with large variations. The highest species diversity (Shannon and Simpson) was found in JM, with 9.71 and 0.9950, respectively, which also showed the lowest values in ACE and Chao1 index, with 2635.51 and 2597.17, respectively (Table 2). The Shannon index value was lowest under the LG, with 9.10, followed by BE (9.20). However, the Simpson index varied considerably and was lowest under PT (0.9902), followed by BE (0.9906). Highest ACE index and Chao1 index value were observed for MB with 3,405.73 and 3,368.12, respectively, and for CB with 3,637.79 and 3,568.24, respectively, most likely because of the high amounts of leaf litter and roots in these sites. In addition, Person’s rank correlation coefficients showed that the Simpson values were significantly positively correlated with the AP (P < 0.05). While, the Shannon value was significantly positively correlated with AK (P < 0.05). Values of ACE and Chao1 were significantly negatively correlated with AK (P < 0.05) (Table 3).

The total number of sequences for each sample in the OTU abundance matrix was randomly sampled at different depths. The rarefaction curve was drawn with the number of sequences extracted at each depth and the corresponding OTU numbers. At a genetic distance of 3%, the rarefaction curve flattened with increasing numbers of measured sequences, indicating that the majority of the sample information was obtained (Supplementary Figure), adequately reflecting the microbial communities in the soil. S2

The relative abundances of bacterial phyla and genera differed among different reforestation pathways (Table 4 and Supplementary Figure S3, S4). Across all samples, we detected 34 bacterial phyla and 1 archaeal phyla, of which 10 phyla had an average abundance greater than 1% (Supplementary Figure S3). The dominant phyla were Proteobacteria (39.98–46.77%), Acidobacteria (15.19–27.32%), Actinobacteria (7.03–16.61%), Verrucomicrobia (4.18–7.97%), Chloroflexi (3.03–7.36%), Bacteroidetes (1.82–4.64%), Gemmatimonadetes (1.93–3.15%), Planctomycetes (1.96–2.48%), Firmicutes (0.21–8.50%), and Nitrospiraea (0.35–3.33%) (Table 4), accounting for 97.12–99.07% of the bacterial sequences from each stands. In addition, Saccharibacteria, Cyanobacteria, Latescibacteria, Elusimicrobia, Armatimonadetes, Chlamydiae, and Parcubacteria were present in most soils, albeit at relatively low abundances (below 1%), and all of the other rarer phyla were identified. Among them, Proteobacteria were the absolute dominant species, followed by Acidobacteria and Actinobacteria (Table 4 and Supplementary Figure S3). The relative abundance of Proteobacteria was highest in PT (46.77%) and lowest in PK (39.98%). The relative abundances of Verrucomicrobia and Planctomycetes in all samples did not significantly different. Actinobacteria were less abundant in coniferous forests [PK (8.26%), LG (7.02%), and PT (8.48%)] than in broadleaf forests [MB (16.40%), JM (15.88%), QM (14.82%), and BE (16.61%)] and conifer-broadleaf forest (14.93%). Firmicutes dominated in PK and were rare in the other stands. While JM had the highest relative abundances of Chloroflexi (7.36%), Gemmatimonadetes (3.15%), and Nitrospirae (3.33%), but the lowest abundances of Acidobacteria (15.19%) and Firmicutes (0.21%) (Table 4 and Supplementary Figure S3).

At the genus level, average relative abundances below 1% were grouped into “others,” and 11 groups were obtained. The relative abundances of several genera differed significantly between the different samples (Supplementary Figure S4). The soil was dominated by Nitrobacter (6.49%), followed by Candidatus-Solibacter (3.72%), Acidothermus (2.86%), Pseudolabrys (2.43%), and Bryobacter (1.98%). Significant differences in relative abundance were not found for the genera Variibacter and Haliangium among samples. The relative abundances of Nitrobacter, Candidatus-Solibacter, Acidothermus, Pseudolabrys, Bryobacter, Rhizomicrobium, H16, and Reyranella showed significant differences between JM and QM. Similarly, the relative abundances of Nitrobacter, Candidatus-Solibacter, Bryobacter, and H16 differed significantly among the JM, PT, and CB. The stand CB had the highest relative abundances of Candidatus-Solibacter, Acidothermus, Bryobacter, and Rhizomicrobium, but the lowest abundance of H16. PT had a significantly higher relative abundance of Nitrobacter compared to the other stands. Significant differences in relative abundance between LG and PK were found for the genera Nitrobacter, Candidatus-Solibacter, Acidothermus, Bryobacter, and H16 (Supplementary Figure S4).

The LEfSe analysis was documented to determine the classified bacterial taxa with significant abundance differences among the different forest stands. As presented in Figure 2, 39 bacterial taxa were showed significantly different with LDA effect size scores were >4. At the phylum level, the biomarkers were affiliated with Proteobacteria, Acidobacteria, Actinobacteria, Chloroflexi, and Bacteroidetes.

We applied heatmap analysis to intuitively display the differences in relative abundances of the top 50 bacterial genera that appeared in all soil samples (Figure 3), which can reflect the differences in soil bacterial community structure between different vegetation types. The resulting heatmap could be divided into five clusters. Both relative abundance and distribution of soil bacteria in different stands were significantly different. In Cluster 1, higher relative abundances of 10 bacterial genera were found in JM, and significant differences were found for the other seven stands. The 19 bacterial genera with higher relative abundances in Cluster 2 were mainly found in QM, BE, and MB. In Cluster 3, the higher relative abundances of the six bacterial genera were mainly found in the PT and CB, and the three bacterial species with higher relative abundances in Cluster 4 were mainly found in PK and were rare in other stands. The three bacterial genera with higher relative abundances in Cluster 5 were mainly found in LG, CB, and MB. This indicates that the compositions and relative abundance of soil bacteria differed significantly between the stands. To further explore these differences, we applied principal components analysis (PCA) to extract the main components. The results of the PCA showed that soil bacterial community from the coniferous forest stand was separated from those of coniferous-broadleaf forest and the broadleaf forest, especially along PC2 (Figure 4).

Canonical correspondence analysis (CCA) was applied to explore the relative abundance of dominant bacteria phyla and genera as a function of soil variables (Figure 5). The first and second axes accounted for 91.2% and 80.5% of the variation, respectively, indicating that soil environmental factors significantly influence the bacterial community structure. At the phylum level (Figure 5A), AP (r = 0.6453) and TK (r = 0.7877) were significantly associated with axis1, which explained 72.4% of the variation. In contrast, soil pH (r = 0.8449) and AK (r = 0.6020) were significantly associated with axis 2. At the genus level (Figure 5B), pH (r = 0.8079) was associated with axis 1, accounting for 55.1% of the variation, while TC(r = 0.6839), TN (r = 0.7548), AN (r = 0.6798), and AK (r = 0.8993) were associated with axis 2. Based on this, soil pH, TC, TN, AN, TK, and AK significantly influenced the bacterial community structure.

We investigated the relationships between the relative abundances of different dominant bacterial phyla and genera and environmental variables (Table 5). In terms of soil chemical properties, AP was negatively correlated with the relative abundances of Proteobacteria, but positively correlated with Bacteroidetes abundance, while Acidobacterias’ relative abundance was negatively correlated with AK (P < 0.05). The relative abundance of Actinobacteria was significantly positively correlated with TC, TN, and AN (P < 0.01), and Chloroflexi relative abundance was significantly positively correlated with pH (P < 0.01). The relative abundances of Firmicutes was positively correlated with TK (P < 0.01).

At the genus level, Acidothermus was significantly negatively correlated with pH, while Nitrobacter was significantly negatively correlated with AP. The abundance of Pseudolabrys was positively correlated with TK, while Candidatus-Solibacter, Bryobacter, and Rhizomicrobium were significantly negatively correlated with AK.

Reforestation pathways greatly affected the soil physicochemical characteristics through various processes, including changes in moisture and temperature, as well as the production of litter and root exudates. In our study, we observed that the broadleaf forests (JM, QM, MB) had higher levels of TC, TN, AN, and AK compared to the coniferous forests (PT, LG, PK) and the conifer-broadleaf forest (CB) (Figure 1), which is consistent with previous research suggesting that broadleaf-Korean pine mixed forest and secondary poplar-birch forest presented significantly higher OC and TN values than spruce fir and larch forest (Hui et al., 2014). The higher TC and TN values in the three broadleaf forests are most likely a result of constant organic matter inputs, rhizodeposition, and the release and recycling of nutrients (Finér et al., 2007; Wang and Cao, 2011). In complete contrast to the TC and TN pattern, the soil C/N ratios in coniferous forests (PT, LG, and PK) and the conifer-broadleaf forest were significantly higher than those in broadleaf forests (QM, BE, and MB). Previous studies have also reported that coniferous forests tend to be highly nitrogen deficient compared with broadleaf forests (Vitousek and Howarth, 1991; Mcgroddy et al., 2004). The soil C/N ratio directly reflects the soil’s nitrogen mineralization capacity, and low soil C/N levels indicate a higher nitrogen mineralization rate, facilitating N uptake by microbes and plants. In turn, high C/N levels are beneficial for the fixation the of soil organic carbon. Litter C/N ratios of coniferous forests are generally higher than those of broadleaf forests, which may contribute to the higher C/N ratio in coniferous forest soils (Yang and Luo, 2011). Several studies have reported that higher C/N ratios and lower nutrient contents in coniferous forests compared with broadleaf forests (Barbier et al., 2008). In our study area, the Liaodong mountainous area, the soil was relatively acidic, with pH value ranging from 4.89 to 5.70, with the lowest levels under QM (pH = 4.89). These differences in pH may be due to different foliage properties and variations in litter quality (Haghdoost et al., 2011). Soil pH, nutrient contents, and numerous biogeochemical processes in forest ecosystems can be influenced by litter chemistry (Mueller et al., 2012; Fanin and Bertrand, 2016). Compared to other broadleaf forests, litter leaf quality in QM was relatively low, with a low nitrogen content, a high C/N ratio, a higher lignin content, and higher lignin/N. Therefore, the decomposition rate of QM litter and the release rates of plant nutrients decreased gradually.

Previous studies have demonstrated the crucial role of soil characteristics in altering the soil microbial communities during vegetation restoration (Thomson et al., 2015). Soil pH and nutrient contents, especially in terms of C, N, and P availability, are paramount factors (Ramirez et al., 2010; Tan et al., 2013) and significantly affect microbial abundance (Christianl et al., 2008; Paulinem and Davide, 2008; Zhong et al., 2010). When soil pH was below 6.5, the pH was supposed to be the primary factor which effected bacterial community diversity, structure and activity, and higher richness and diversity values were found at a pH close to neutral (Preem et al., 2012; Ramirez et al., 2012; Bergkemper et al., 2016). In contrast, Lauber et al. (2009) has proposed that there is a significantly negative relationship between soil pH and bacterial diversity. However, in my study, soil pH was not significantly correlated with any bacterial diversity index. However, we observed a relatively small pH rangs (4.89 to 5.70), making it different to find any correlation with bacterial diversity. Previous studies have suggested that bacterial alpha diversity is largely affected by SOC and TN (Siles and Margesin, 2016; Ren et al., 2017). However, in our study, TC and TN levels were not significantly correlated with bacterial diversity indexes, which was consistent with the observations obtained by Lin et al. (2014).

In most terrestrial ecosystems, the soil C/N ratio greatly influences soil microbial communities (Fierer and Jackson, 2006) and is negatively correlated with the soil bacterial Shannon index (Zhou et al., 2017). However, in our study, C/N ratio showed no correlation with bacterial diversity indices, most likely because of the relatively small C/N range (10.64–13.81).

However, the Simpson index was significantly positively correlated with the AP, which is in accordance with a previous study suggesting that P increases bacterial diversity in pasture soil (Tan et al., 2013). Similarly, soil P is the second most significant driver of bacterial diversity in soil (Siciliano et al., 2014), although some authors have stated that bacterial diversity was not affected by P levels (Liu et al., 2018). In our study, the Shannon value was significantly positively correlated with AK (P < 0.05) (Table 3) and increased ACE index and Chao1 index values were observed in CB and MB, reflecting the large amounts of leaf litter and roots. Both of these α-diversity indicators were significantly negatively correlated with AK (Table 3), which leads us to infer that AP and AK availability may be the limiting factor affecting bacterial diversity in this area.

Similar to bacterial community diversity, the relative abundances of dominant bacterial phyla were influenced by different revegetation types. With regard to bacterial compositions, the different revegetation forest stands harbored distinct bacterial communities. In our study, Proteobacteria, Acidobacteria, and Actinobacteria were the most abundant soil bacterial phyla in the Baishilazi Nature Reserve, which is in agreement with a past study in northern soils (Sun et al., 2014). In the light of bacterial classification, Proteobacteria, Acidobacteria, Actinobacteria, Firmicutes, Gemmatimonadetes, Nitrospirae, and Chloroflexi significantly differed between the different forest stands (Table 4). Acidobacteria was the dominant taxa in LG, while Actinobacteria was the dominant taxa in PK; Proteobacteria dominated the soil under PT. In constrast, Nitrospirae and Chloroflexi were dominant in JM. In contrast to our study, previous studies reported that coniferous forests were dominated by Proteobacteria and Gemmatimonadetes (Hui et al., 2014).

Proteobacteria was the most abundant phylum in our research, which is similar to the result of previous studies (Li et al., 2014a; Miyashita, 2015). However, according to Liu et al. (2013), Actinobacteria play a dominant role in the bacterial community of the mixed forest soil of the Dinghushan Mountain. Generally, members of the Proteobacteria and Acidobacteria are ubiquitous in almost all soil types (Zhang and Xu, 2008). In our study, despite different revegetation types, including broadleaf forests, coniferous forests, and conifer-broadleaf forest, the composition of bacterial groups was similar for all phyla detected. Previous studies have suggested that the Proteobacteria can be used as indicators of the prevalent nutrient status due to their different natures (Hartman et al., 2008). In our study, the relative abundance of Proteobacteria was highest in PT and largely affected by AP, but not by measured C, although a high abundance of Proteobacteria has previously been associated with a high C availability (Fazi et al., 2005; Fierer et al., 2007). This may be related to the composition of Proteobacteria in the different forest stands, and our results indicated that in the study region, soil TC had no significant impact on the bacterial community composition.

Recent studies have shown that Acidobacteria are generally oligotrophic (Kielak et al., 2009; López-Lozano et al., 2013), and Acidobacteria play a key role in biogeochemical nutrient cycling (Tringe et al., 2005). Apparently, the level of nitrogen input is negatively correlated with the relative abundance of Acidobacteria (Fierer et al., 2012), which is line with the results of our study. Acidobacteria belonged to the acidophilic bacteria, and Acidobacteria abundance is significantly higher under acidic conditions (Dion, 2008; Lauber et al., 2009; Bardhan et al., 2012). However, in our study, the relative abundance of the Acidobacteria was not correlated with soil pH (Table 5) but was significantly negatively correlated with AP. Soil properties, especially AP and AK, were significantly correlated with the bacterial, even at the very small geographical scale in our study, and this finding is in line with the previous research (Cai et al., 2018).

In our study, except for the Proteobacteria and Acidobacteria, the significant difference in the relative abundance of Actinobacteria was observed in coniferous forests (PK, LG, and PT) and was lower than the abundance in broadleaved forests (JM, QM, and BE). This finding is in agreement with a previous which reported that Actinobacteria are among the most important organic matter decomposers in soil (Kopecky et al., 2011). In our study, the differences in relative Actinobacteria abundance between the different reforestation pathways were related to changes in TC, TN, and AN (Figure 5 and Table 5), which is in agreement with a previous study stating that Actinobacteria are significantly positively correlated with the soil organic matter accumulation, during secondary forest succession (Guo et al., 2018). In addition, PK showed a higher relative abundance of Firmicutes compared to the other forest stands, confirming the high resistance of these bacteria to unfavorable conditions (Hartmann et al., 2014). Other specific taxa, particularly Chloroflexi and Nitrospirae, were significantly correlated with the soil pH (Table 5).

The clear differentiation of bacterial communities and the clustering in the heatmap (Figure 3) and in PCA (Figure 4) plots, suggests that broadleaf forests and coniferous forests harbored different bacterial communities, indicating that the soil bacterial community composition largely depends on soil physicochemical characteristic. Soil physicochemical variables therefore play the decisive role in altering bacterial communities during vegetation restoration, which is consistent with previous studies (Kuramae et al., 2010; Li et al., 2014b).