A comparison of microbial composition under three tree ecosystems using the stochastic process and network complexity approaches

The study area is located in the Liupanshan National Nature Reserve (106°09′-106°30′E, 35°15′-35°41′N), which is in the successional zone of the warm-temperate semi-humid and semi-arid region in northern China, with high temperatures and rain in summer and dry and cold in winter. The average elevation temperature is 17.4°C, and the average precipitation is 676 mm. The coniferous forests in the region are mainly L. principis-rupprechtii, and the broad-leaved forests are mainly Quercus wutaishanica, Betula platyphylla, Betula albo-sinensis, and Populus davidiana. In the nature reserve, nine replicate plots of 20 × 20 m (27 in total) were designated as coniferous forest (CF) (9 plots), coniferous and broad-leaved mixed forest (MF) (9 plots) and broad-leaved forest (BF) (9 plots) at similar elevations, and the distance between adjacent plots was greater than 200 m (Figure 1). According to the sampling method of Pan et al. (2021), soil samples from each plot were collected by 5-point sampling method after removing litter from the surface. Plant roots were avoided during sampling and soil samples were collected at a depth of 10–20 cm. Five soil samples were mixed and defined as the soil sample of one plot, which was then sieved with 2 mm and divided into two parts. One part was stored in dry ice to extract soil nucleic acids, and the other was used for soil physicochemical characterization (Lin et al., 2021).

Soil pH was measured by pH meter (Mettler S220, Mettler Toledo Solutions, Greifensee, Switzerland) (Wang et al., 2019). And soil water content (SWC) was determined using the weighing method before and after drying at 105°C (Bao, 2000). After the 27 collected soil samples were air-dried, the soil total carbon (TC), total nitrogen (TN), and total phosphorus (TP) were determined using atomic absorption spectrometry (iCE 3,500, Thermo Fisher Scientific, Waltham, MA, United States) (Bettinelli et al., 2000).

The total genomic DNA from 27 soil samples was extracted by CTAB method according to the manufacturer’s instructions. The quality and quantity of the extracted DNA were measured using a NanoDrop™ 2000 spectrophotometer (Thermo Fisher, United States) and 1% agarose gel electrophoresis, respectively (Lian and Xing, 2017). The soil microbial community was examined using Illumina NovaSeq 6,000 sequencing analysis (Illumina, San Diego, CA, United States). The Illumina NovaSeq 6,000 sequencing libraries for bacteria were prepared via PCR amplification of the V3-V4 hypervariable regions of the bacterial 16S rRNA gene using the primers 338F (5’-ACTCCTACGGGAGGCAGCA-3′) and 806R (5’-GGACTACNNGGGTATCTAAT-3′) with the GeneAmp® 9,700 PCR thermocycler (Applied Biosystems, CA, United States) (Caporaso et al., 2012). The Illumina MiSeq sequencing libraries for fungi were prepared via PCR amplification of the internal transcribed spacer (ITS) region of fungi using the primers, ITS5 (5’-GGAAGTAAAAGTCGTAACAAGG-3′) and ITS2 (5’-GCTGCGTTCTTCATCGATGC-3′) (Schmidt et al., 2012). The thermal cycling conditions consisted of an initial denaturation at 95°C for 3 min followed by 27 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 45 s, followed by a final extension at 72°C for 10 min (Haas et al., 2011). The purified amplicons were pooled in equimolar concentrations and paired-end sequenced on Illumina platform according to the standard protocols of Novogene Bioinformatics Technology Co. Ltd. (Beijing, China).

Paired-end (PE) library were constructed using a NEXTflex Rapid DNA-SEQ Kit (Bioscience, South San Francisco, CA, United States), and an Illumina NovaSeq 6,000 platform was used for sequencing. Trimmomatic software (Illumina) was used for quality control of Illumina MiSeq sequencing original sequences. FLASH (1.2.11)1 software was used for stitching (Magoč and Salzberg, 2011). UPARSE 7.1 software was used for amplicon sequence variants and operational taxonomic units (OTU) clustering analysis (similarity 97%), and UCHIME software was used to remove chimeras (Caporaso et al., 2012). Each sequence was annotated for species classification by the ribosomal database project classifier and compared with the Silva database (SSU128) at a confidence threshold of 0.7 (Quast et al., 2013).

The QIIME program (1.9) was used to calculate the alpha diversity indices (Shannon, and ACE) (Kuczynski et al., 2012). The “linkET”2 package in R software (Version 4.1.0) was used to calculate and visualize the correlation between soil physicochemical properties and soil bacterial alpha diversity. A non-metric multidimensional scale (NMDS) based on the Bray-Curtis distance matrix and ANOSIM test with 9,999 permutations were used to illustrate beta diversity (Ziegler et al., 2017). Differences in community composition between sample groups were analyzed using the “vegan” package, and the “Hmisc” and “minpack.lm” package in R software were used to construct the neutral community model (NCM) of bacterial and fungal communities (Chen et al., 2019). The “NST” package in R software was calculated the βNTI and RC values for bacterial and fungal community. The soil physicochemical properties and soil bacterial and fungal communities (Spearman correlation) were calculated based on the Bray–Curtis dissimilarity matrices used “vegan” package (Wang et al., 2017). The “ggalluvial” package in was used to describe changes in bacteria and fungi (top 5 phyla with the highest abundance), and the “psych” package was used to calculate the correlation between top 5 phyla and soil properties, which was then visualized using the “ggplot2” package.

In order to better understand the connectivity and complexity of bacterial and fungal communities in CF, MF, and BF sample plots, we used “psych” package to calculate the Spearman correlation coefficients among OTUs in all samples, and the detection rate was greater than 77.7%. At the same time, the false discovery rate was used to correct the value of p multiplicity test. The rank correlation coefficient r > |0.9| and p < 0.001 were determined (Berry and Widder, 2014). Then, Cytoscape (3.7.1) was used for network visualization, and the number of OTUs at the phylum level was statistically analyzed. The top 6 modules with the highest OTUs enrichment in each network were selected for coloring. These modules were extracted by Cytoscape, colored with dominant phylum, and visualized using “ggplot2” package. The “psych” package was used to calculate the correlation between the OTU with the top 20 degrees in bacteria and fungi and soil properties, which was then visualized using the “ggplot2” package.

SPSS version 25.0 (SPSS Inc., Chicago, IL, United States) was used for one-way analysis of variance of soil physicochemical properties data, and Duncan’s multiple range test was used to identify the significant differences between means at a 5% significance level. All data are presented as the mean ± standard error (n = 9; Pan et al., 2021).

With the succession from CF plot to BF plot, soil pH, SWC, TC and TP contents increased, and BF plot was significantly higher than CF plot (p < 0.05). It was worth noting that soil TN content in the MF plot was the highest and significantly higher than that in the CF plot. Whereas, soil C:N in the BF plot was higher than that in the CF and MF plots, respectively. In addition, soil N:P and C:P changed in the MF plot were similar to those in TN (Figure 2A).

Alpha diversity analysis showed that the number of bacterial OTUs and bacterial richness in the three sample plots were not significantly different, but the soil bacterial diversity index was higher in the MF sample plot and significantly higher than in the CF and BF sample plots. Interestingly, there was no significant difference in soil fungal community diversity among the three plots, while the number of fungal OTU and fungal richness were consistent, that is, with the succession from CF plot to BF plot, the number of fungal OTU and fungal richness increased, and BF plot was significantly higher than CF plot (p < 0.05; Figure 2B).

Further, the Mantel test showed that soil fungal richness in CF plot was significantly correlated with C:N (p < 0.05). The number of bacterial OTU was strongly correlated with pH, bacterial diversity was strongly correlated with C:N, and bacterial richness was strongly correlated with TC and C:P. With the succession from CF to BF plot, there was a correlation between the alpha diversity index of soil bacteria and fungi and the soil physicochemical properties (p < 0.05), such as bacterial OTUs and SWC, bacterial diversity and TP, bacterial richness and TN, and fungal richness and N:P index (Figure 2C).

NMDS analysis showed that the composition of the soil bacterial (R = 0.6491, p = 0.0001) and fungal (R = 0.4966, p = 0.0001) communities were different among the three plots (Figure 3). NCM analysis showed that stochastic processes gradually dominated the soil bacterial community composition with succession from CF to BF plot (R² = 0.785–0.831; Supplementary Figure S1). Further analysis by βNTI and null model showed that with the succession from CF to BF plot, the ecological processes of soil bacterial communities changed from deterministic processes (69.4% in CF plot) to stochastic processes (77.8% in BF plot). In CF plot, homogeneous selection (55.6%) was dominant; while in BF plot, drift (55.6%) was dominant. The ecological process of soil fungal community was similar to that of bacteria. However, homogeneous selection (86.1%) was dominant in CF plot, while heterogeneous selection (52.8%) and drift (47.2%) were dominant in BF plot (Figure 3).

In order to further explore the relationship between Bray–Curtis dissimilarity of soil bacteria and fungi and environmental factors, we found that Bray–Curtis dissimilarity of soil bacteria was positively correlated with soil pH, SWC and TN (Figure 4A). Fungal Bray-Curits dissimilarity was significantly positively correlated with soil SWC, TC and TN (p < 0.05; Figure 4B). The Bray–Curtis dissimilarity of soil bacteria and fungi were significantly correlated with environmental factors (p < 0.05).

Of the 27 soil samples, 84.57% OTUs of bacteria and 36.40% OTUs of fungi were annotated at the phylum level. With the succession from CF to BL plot, the relative abundance of Proteobacteria (29.2–27.6%), Actinobacteriota (7.1–5.3%) and Bacteroidota (4.1–4.0%) decreased, while the relative abundance of Acidobacteriota (10.7–11.2%) and Verrucomicrobiota (4.6–5.0%) increased (Figure 5A). The relative abundance of Ascomycota (33.1–32.3%), Mortierellomycota (22.2–14.9%), Rozellomycota (3.6–1.1%) and Chytridiomycota (1.3–1.1%) decreased, while the relative abundance of Basidiomycota (14.8–21.2%) increased (Figure 5B).

The correlation between the relative abundance of dominant phyla and the soil physicochemical properties was further analyzed. Actinobacteriota were significantly negatively correlated with soil pH, SWC, TC and TN (p < 0.05), while Proteobacteria and Verumobacteria were significantly negatively correlated with TP and TN, respectively (Figure 5C). Rozellomycota was significantly negatively correlated with soil pH, SWC and TC, while Mortieromycota was significantly negatively correlated with TP (p < 0.05; Figure 5D).

The top 20 bacterial genera in relative abundance of three sample plots were MND1, Candidatus Udaeobacter, Sphingomonas, RB41, Brevundimonas, Haliangium, Escherichia Shigella, Bradyrhizobium, IS44, Rhodoplanes, Ellin6055, Candidatus Solibacter, Ellin6067, Dongia, Bryobacter, Agathobacter, Gaiella, Terrimonas, Pedomicrobium and Subgroup 10 (Supplementary Figure S2A). And MND1, RB41, Candidatus Solibacter, Ellin6067, and Subgroup 10 were the differential genera of top 20 bacterial genera in relative abundance in CF and BF sample plots (Supplementary Figure S2B). While the top 20 fungal genera in relative abundance of three sample plots were Mortierella, Fusarium, Inocybe, Tuber, Pseudogymnoascus, Penicillium, Cladosporium, Hymenogaster, Malassezia, Pichia, Trichophaea, Trichocladium, Aspergillus, Solicoccozyma, Chaetomium, Protoglossum, Archaeorhizomyces, Phallus, Geminibasidium and Dactylonectria (Supplementary Figure S3A). However, there were no differential genera of top 20 fungal genera in relative abundance in CF and BF sample plots (Supplementary Figure S3B).

The co-occurrence network showed that the network connectivity and complexity differed between soil bacteria and fungi. In the bacterial network, the network relationship in the MF plot was more complex (1,818 nodes and 3,516 edges), and the BF plot (1,798 nodes and 2,831 edges) had higher network connectivity and complexity than the CF plot (1,624 nodes and 2,676 edges). The number of soil bacterial network modules in CF plot (105) was also higher than that in MF (87) and BF (88) plots. In addition, the positive edges, negative edges and average degree of soil bacterial network in MF (2,573 positive edges, 943 negative edges) were higher than those in CF (1,950 positive edges, 729 negative edges) and BF (2,122 positive edges, 709 negative edges) plots. It was worth noting that the BF plot had the highest positive/negative edge (2.99). The average path distance of BF plot (10.67) was higher than that of CF (10.24) and MF (8.79), and the harmonic geodesic distance of BF plot (9.26) was also higher than that of CF (8.46) and MF (7.61; Figure 6; Supplementary Table S1).

In the soil fungal network, the connectivity and complexity of the network gradually increased with the succession of CF plot to BF plot (the number of nodes from 47 to 149, and the number of edges from 100 to 125). The number of soil fungal network modules in BF (65) plot was also higher than that in CF (37) and MF (12) plots (Figure 7; Supplementary Table S1).

In the network analysis, the OTUs in the main modules of soil bacteria in the three sample plots were dominated in Acidobacteriota, Proteobacteria, Actinobacteriota, and Bacteroidota; in fungi, Ascomycota, Basidiomycota and Mortierellomycota were predominant. The top 20 OTUs of each network in terms of degree were selected, and correlation analysis with soil physicochemical properties showed that the magnitude of the bacterial network degree was closely related to physicochemical properties. Among them, soil bacteria in CF showed a negative correlation with pH and a positive correlation with TC; soil bacteria in MF showed a negative correlation with TC, TN, and C:P ratio; and in BF, soil bacteria showed a positive correlation with TC, C:N, N:P and C:P ratio (p < 0.05; Supplementary Figure S3; Figure 8).

In this study, the results showed significant soil pH, SWC, TC, TN, and TP differences among CF, MF, and BF (p < 0.05; Figure 2A). Differences in plant species have been reported to lead to changes in soil properties caused by species differences affecting root exudates and litter, resulting in an unbalanced input of soil nutrients (Li C. Y. et al., 2021). In general, the litter of BF contributes significantly to soil C, N, and P accumulation (Pang et al., 2020). In this study, the soil TC, TN, and TP contents in BF plot were higher than those of MF and CF. Previous study pointed out that the difference in litter quality, the leaves of broad-leaved forest were larger and decomposed faster than those of coniferous forests, which was conducive to the accumulation of nutrients in the soil, whereas coniferous forests were the opposite (Mueller et al., 2015). Soil pH and TC content increased in the succession from coniferous forest to broad-leaved forest (Qu et al., 2020), and similar conclusions were obtained in this study. Interestingly, the MF maintained high N:P ratio. The existing studies suggested that the leaves of trees, shrubs and herbaceous plants in mixed forests maintain high N:P ratios (Pang et al., 2020). With the input of plant leaf litter to soil nutrients, the reason why N:P was higher in MF plot than CF and BF in this study was speculated. It was noteworthy that the N:P ratios of the three sample plots were all less than 10, indicating that soil in Liupanshan Mountain was obviously limited by nitrogen (Elisabeth and Brent, 2013).

Vegetation type can indirectly change the diversity and richness of soil bacterial and fungal communities by changing the soil physicochemical properties (Yao et al., 2017; Deng et al., 2019). In this study, there were differences in soil bacterial community diversity among the three sample plots (Figure 2B). A recent study indicated that the coniferous and broad-leaved mixed forest had higher soil bacterial alpha diversity, which was due to the higher soil TN (Li et al., 2014). In addition, the richness of soil fungal community increased with the succession from CF to BF plot (Figure 2B). It has been reported that coniferous and broad-leaved tree species mixing can increase the richness of soil fungal community and subsequently affect the structure of soil fungal community, which further demonstrated that differences in soil physicochemical properties (temperature, soil moisture, litter, TN, C:N and N:P ratio) were important factors affecting the structure of soil microbial community (Dong et al., 2021).

According to Mantel test, soil fungal community richness was significantly correlated with soil C:N in CF plot (p < 0.05; Figure 2C). Because coniferous forest had more difficult to degrade lignin and coarse fiber, these substances were difficult to be broken down by microorganisms. Since soil fungi participate in the decomposition of recalcitran organic matter (lignin), fungi that were more dependent on soil nutrients have a strong correlation with C:N (de Boer et al., 2005; Šnajdr et al., 2008). This is consistent with earlier findings and provides new evidence revealing the effects of tree species differences on soil properties and microbial communities in temperate forests (Gao et al., 2019; Sheng et al., 2019). An increasing number of studies have recorded the response of soil microbial community diversity to environmental changes (de Carvalho et al., 2016; Tripathi et al., 2016). Among them, soil bacterial diversity and community composition were observed to be strongly correlated with soil pH (Rousk et al., 2010).

Deterministic and stochastic ecological processes represent complementary and continuous ecological processes of assembly patterns in microbial communities (Zhou and Ning, 2017). In our study, the driving effect of stochastic processes on soil bacterial and fungal communities gradually increased as the sample plots changed from CF to MF (Supplement Figure S1). Therefore, we infer that there is a unified continuum between deterministic and stochastic processes in the microbial community in our study area, which is consistent with the results of previous studies (Chase et al., 2011; Tian et al., 2017). Myers et al. (2013) indicated that the assembly patterns of soil microbial communities in temperate forests reflected stronger environmental correlations, whereas plant species differences played a key role. With succession from CF to BF, the stochastic process of the soil bacterial community gradually dominated. The study further showed that the process of bacterial community assembly patterns also changes with the long-term change of vegetation cover (Goss-Souza et al., 2017). In the early stages of succession, the assembly patterns of soil bacterial communities were dominated by stochastic processes, whereas in the late stages of succession, deterministic processes gradually dominated (Ferrenberg et al., 2013; Dini-Andreote et al., 2015). In this study, L. principis-rupprechtii was planted early and for a long time, confirming the above conclusions. In addition, the succession from CF to BF was driven but not dominated by stochastic processes in the soil fungal community (Figure 3). One explanation for the differences in the assembly patterns of bacterial and fungal communities is that fungal communities are initially strongly dominated by deterministic processes; however, this dominance gradually diminishes as succession proceeds along the temporal timescale (Tian et al., 2017). Nonetheless, soil bacterial communities have more remarkable metabolic functional plasticity (Massana and Logares, 2013), which may be because the deviation level of soil fungal communities is lower than that of soil bacterial communities. This also explains our first hypothesis that stochastic and deterministic processes influence community assembly patterns of soil bacteria and fungi in temperate forests.

Environmental heterogeneity was the main factor affecting the geographical distribution of the microbial communities (Zhou et al., 2013; Liu et al., 2019). For example, soil salinity, pH, carbon and nitrogen content, and the C:N ratio may all be determinants of the soil microbial community (Griffiths et al., 2011; Xiong et al., 2014). In this study, the soil bacterial and fungal communities were affected by environmental factors. Soil pH and SWC content were strongly correlated with the assembly patterns of the soil bacterial communities (Figure 4). Previous studies have shown that the degree to which stochastic and deterministic ecological processes in the soil bacterial community are affected by soil pH (Tripathi et al., 2018), and SWC also explains the changes in soil bacterial communities (Delgado-Baquerizo et al., 2018). Soil SWC and TN content were significantly correlated with fungal communities (p < 0.05), and previous studies have shown that changes in bacterial communities largely depended on changes in soil pH (Rousk et al., 2010). In contrast, the main influences on fungal communities tended to be soil carbon, moisture, and plant species and were not sensitive to changes in soil pH (Barberán et al., 2015; Prober et al., 2015). In our study, the soil TN content of CF was significantly different from that of MF and BF, which highlighted the critical role of soil physicochemical properties in fungal community assembly.

Soil bacterial and fungal communities with different ecological functions play important roles in forest ecosystems. In our study, Actinobacteriota had a strong correlation with pH, SWC, TC and TN (Figure 5). It has been shown that Actinobacteriota, as k-strategists oligotrophic bacteria, depends on the availability of water and nutrients (Li et al., 2022), and can play an important role in microbial co-occurrence networks (Chen et al., 2021). Our study further confirmed the above conclusions. In addition, a study on forest fungal diversity concluded that changes in total soil carbon and nitrogen content during forest succession were negatively correlated with the relative abundance of Rozellomycota (Bai et al., 2019), which explained the decrease in the relative abundance of Rozellomycota during the succession from CF to BF sample plots in this study and the significant correlation with soil TC and C:P. It has been hypothesized that an increase in the relative abundance of Mortierellomycota decreases soil microbial diversity (Li X. et al., 2021), and our study found that the relative abundance of Mortierellomycota decreased with the succession from CF to BF sample plots and was significantly correlated with soil TP content (p < 0.05). These results further proved that the change of soil available nutrients during forest succession was determined by the interaction between microorganisms and litter inputs of different quality.

Network analysis to study microbial co-occurrence patterns can better reveal interactions between ecological niches and microbial taxa (Barberán et al., 2012). We constructed bacterial and fungal co-occurrence networks for each of the three sample plots and found that the microbial co-occurrence networks did not occur randomly. However, communities or OTUs formed modular groupings of differently associated species (Barberán et al., 2012; Sul et al., 2013). In network analysis, a higher degree of modularity and denser network recorded a highly complex microbial community (Olesen et al., 2007). In this study, the number of soil bacterial network modules was the largest in CF, and the number of soil fungal network modules was the largest in BF (Figures 6, 7), which has been confirmed in studies related to soil bacterial and fungal community co-occurrence networks during forest ecosystem restoration (Sun et al., 2017); that is, the bacterial co-occurrence networks became more complex as succession proceeded (Avera et al., 2015; Xue et al., 2017). In addition, the number of nodes and edges in the network further characterizes the connectivity of the network. The nodes and edges of soil bacteria in CF were lower than those in MF and BF, and the positive/negative edge ratios were the lowest in this study (Supplementary Table S1). This implies that in the co-occurrence network, microorganisms compete for ecological niches favorable to their growth needs through competition and complementarity (Faust and Raes, 2012). During grassland restoration, microorganisms avoid grazing intrusion through ecological niche competition to maintain the soil ecosystem function (Yao et al., 2018; Cong et al., 2021). In another study, with seasonal changes, plant rhizosphere bacterial communities responded to utilizing and acquiring soil nutrients during the plant growth stages through a complementary ecological niche strategy (Shi et al., 2016). The average path distance in the network highlighted the shortest path between two nodes. From an ecological point of view, the number of modules and complexity of paths characterizes the cluster species as performing more functions and roles in the network (Wang et al., 2016). In this study, the fungal communities in the CF differed from those in the MF and BF and had more complex network relationships, indicating that the fungal communities in the CF were more competitive. In general, abundant resources favored taxa selection and led to enhanced interactions among OTUs in the network, whereas soil nutrient limitation intensified resource grabbing among taxa and formed closer interrelationships (Pan et al., 2021).

We further analyzed the correlation between the top 20° of the major taxa in each network and the physicochemical properties, and found that with the succession from CF to BF plots, the correlation between the major taxa in each network and the physicochemical properties of soil was different. For example, bacterial taxa in the CF plot were almost negatively correlated with pH, but strongly positively correlated with TC (p < 0.05; Figure 8). Members such as Gaiella and RB41 are known to be associated with inputs to soil pH and carbon inputs (Lin et al., 2019; Wang et al., 2022). This further confirmed the influence of refractory carbon input in coniferous forests on soil microbial interactions. In the MF plot, members of Vicinamibacterales were mostly negatively correlated with soil TC and TN. In view of the potential role of Vicinamibacterales in biological antagonism and ecological restoration (Chun et al., 2021; Xu et al., 2022); it is speculated that they can be used as indicator taxa of mixed temperate forest and regulated by soil carbon and nitrogen. Major taxa such as Afipia, Rhodopseudomonas and Mesorhizobium actively participated in the soil carbon, nitrogen and phosphorus cycling, and the significant correlation with physicochemical properties in the BF plot further validated the results of previous studies (LaSarre et al., 2017; Duan et al., 2020; Luo et al., 2020). In addition, the dependence of soil fungal communities on carbon resources in BF plot further indicated that the interactions of major taxa could further stabilize the microbial community structure. This also confirmed our second hypothesis that succession from CF to MF and BF plots changed the network relationship of soil bacterial and fungal communities.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/↗, PRJNA867578, PRJNA868111.