Similarities and differences in the microbial structure of surface soils of different vegetation types

Site information

The sites were situated in Longzhou County, Chongzuo City, Guangxi Zhuang Autonomous Region, China (106°33′11″–107°12′43″E, 22°8′54″–22°44′42″N). It is characterized by a southern subtropical monsoon climate, characterized by high temperatures, abundant rainfall, and ample sunshine throughout the year. The region experiences a consistent pattern of hot and dry seasons, with approximately 350 frost-free days annually and a frost period lasting for 13 days. Four representative vegetation types were selected for the experiment, namely a woodland with the dominant tree species Horsfieldia hainanensis (HH), a woodland with the dominant tree species Drypetes perreticulata (DP), a Zea mays farmland (ZM) and a Citrus reticulata farmland (CR). Both farmland areas were natural forests that were deforested to establish agricultural land. The Zea mays farmland has been under cultivation for 20 years, while the Citrus reticulata farmland has been cultivated for 10 years. Both farmland areas have been regularly fertilized as part of their daily management.

Soil sampling

In December 2022, soil samples of HH, DP, ZM and CR were taken from 0 to 10 cm. According to the S-type sampling principle, eight soil cores of 0-10 cm were randomly selected in each plot and mixed to obtain a total of 16 samples. Each composite soil sample was carefully collected and placed in a plastic bag. The bag was labeled to ensure proper identification of the sample. Subsequently, the samples were transported to the laboratory in an ice box to maintain their integrity. In the lab, stones and plant residues such as roots and litter were removed from the soil. The soil material was then passed through two mm sieve to remove any coarse particles. Following this, the sieved soil was immediately transferred into 2 ml centrifuge tubes and frozen at a temperature of −80 °C. This freezing process was carried out to facilitate the extraction of microorganisms from the soil at a later stage.

Soil chemical properties determination

10g of soil sample was weighed and placed in a 50 mL conical flask, and double distilled water was added according to the principle of soil/water ratio of 2.5:1. After 2 min of high-speed shaking, the soil pH value was measured for half an hour by using a pH meter (Lei-ci PHS-3C). Dried soil samples of 10–13 mg were weighed and sealed in a tin container and then the total carbon and nitrogen contents of the soils were determined by the elemental analyzer (Elementar Vario EL III Germany). Using a molybdenum antimony anti-colorimetric method and HClO₄-H₂SO₄ heating digestion, the total phosphorus content of the soil was determined. 1 g of the sample was to be weighed in a cleaned and dried digestion tube. eight mL of concentrated sulfuric acid (H₂SO₄) were then required to be mixed, and allowed to soak for an overnight period. Finally, 10 drops of HClO₄ were supposed to be put in. The digestion tube might then be heated on the digestion machine until the boiling liquid is clarified. The sample was kept in a 100 mL volume bottle after digestion and cooling. The molybdenum-antimony resistance colorimetric method was used to perform colorimetric analysis on five mL of filtrate using a spectrophotometer (P4 UV-Visible China) and a 50 mL volumetric bottle (Liu et al., 2022b).

DNA extraction and high-throughput sequencing

The OMEGA Soil DNA Kit (M5635-02) (Omega Bio-Tek, Norcross, GA, United States) was used to extract DNA from all the samples. A NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) was used to measure the quantity and quality of DNA and agarose gel electrophoresis for extracted DNA (agarose concentration of 1.2%). With the assistance of the primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′), the 16S V3_V4 region of the soil bacterium was amplified (Claesson et al., 2009). Additionally, the fungal ITS region was amplified using primers ITS5 (5′-GGAAGTAAAAGTCGTAACAAGG-3′) and ITS2 (5′-GCTGCGTTCTTCATCGATGC-3′) (White et al., 1990). The PCR system had a total volume of 25 µL, consisting of the following components: 5 µL of 5 × reaction buffer, 5 µL of 5 × GC buffer, 2 µL of dNTP (2.5 mM), 1 µL of forward primer (10 µM), 1 µL of reverse primer (10 µM), 2 µL of DNA template, 8.75 µL of ddH2O, and 0.25 µL of Q5 DNA polymerase. The amplification process involved an initial denaturation at 98 °C for 2 min, followed by cycling at 98 °C for 15 s, 55 °C for 30 s, and 72 °C for 30 s for a total of 25 cycles. A final extension was performed at 72 °C for 5 min, with a 10 °C hold for the 25 cycles. To purify the PCR amplicons, Vazyme VAHTSTM DNA Clean Beads from Vazyme in Nanjing, China were utilized. The Quant-iT PicoGreen dsDNA Assay Kit from Invitrogen in Carlsbad, CA, United States was employed for quantification. 250 pair-end sequencing was carried out using the Illlumina NovaSeq platform and NovaSeq 6000 SP Reagent Kit (500 cycles). Amplicons were pooled in equal amounts following the individual quantification step. The above operations were completed in Shanghai Personal Biotechnology Co., Ltd.

Sequence analysis

Microbiome bioinformatics analysis was carried out using QIIME 2 (2019.4), with a slight modification based on the methodology described by Bolyen et al. (2019). The initial step involved the removal of primers using the cutadapt plugin, following the demultiplexed of raw sequence data using the demux plugin (Martin, 2011). Subsequently, the DADA2 plugin was employed to perform quality filter, denoising, merging, and removal of chimera from the obtained sequences (Callahan et al., 2016). To construct a phylogenetic tree, non-singleton amplicon sequence variants (ASVs) were utilized in conjunction with the fasttree2 and mafft tools, which were implemented as part of the pipeline (Katoh et al., 2002; Price, Dehal & Arkin, 2009). Alpha-diversity metrics including Chao1 (Chao, 1984), Shannon (1948), and Pielou’s evenness (Pielou, 1966), were estimated by employing the diversity plugin. In addition, beta diversity metrics (Bray-Curtis dissimilarity) were calculated. The feature-classifier plugin was utilized to assign taxonomy to the ASVs using the naive Bayes classifier and two reference databases: SILVA Release 132 Database for bacteria and UNITE Release 8.0 Database for fungi, which were selected based on the studies conducted by Kõljalg et al. (2013) and Bokulich et al. (2018), respectively.

Data analytics

The main tools used to analyze sequence data were QIIME2 (2019.4) and R packages (v3.2.0; R Core Team, 2015) (New Zealand). Using the ASV table in QIIME2 (2019.4), alpha diversity indices at the ASV level, including the Chao1 richness index, Shannon index, and Pielou’s evenness, were calculated and displayed as box plots. To visualize the differences in alpha diversity between the various specimen groups, the data in the table above were plotted as box plots using QIIME2 (2019.4) and the ggplot2 package for the R package (v3.2.0; R Core Team, 2015). The significance of the differences could be confirmed using the Kruskal–Wallis rank sum test and the dunn’test as a post hoc test (Wickham, 2016). Based on the occurrence of ASVs across samples/groups regardless of their relative abundance, a Venn diagram was created to visualize the shared and unique ASVs among samples or groups using the R package “VennDiagram” (Zaura et al., 2009). PERMANOVA (Permutational Multivariate Analysis of Variance) was used to evaluate the significance of the differences in microbiota structure between groups (Anderson & Willis, 2003). Using abundance information from the top 20 orders of average abundance, the heatmap was produced using R’s pheatmap package. A traditional multidimensional scaling (cMDScale) analysis technique is principal coordinates analysis(PCoA) (Ramette, 2007). This was accomplished by maximising the distance relationships between the original samples and expanding the sample distance matrix in low-dimensional space following projection. The number of permutation tests used in the “permanova” analysis of variance between groups was set to 999. The analysis of variance between groups was performed using the scikit-bio package in Python. After removing singletons from the feature list, the QIIME2 (2019.4) “qiime taxa barplot” was used to visualize the compositional distribution of each sample at the taxonomic levels of phylum and genus. In order to better understand the relationship between soil nutrient content and soil microbial community composition (phylum and genus level), we performed correlation heat map analysis of them based on Spearman rank correlation coefficient. This analysis was performed by the genescloud tools, A free online platform for data analysis (https://www.genescloud.cn) (Liu et al., 2022a). The differences between the chemical characteristics of the soils of the various vegetation types were examined using a one-way ANOVA test. Waller-Duncan post-hoc multiple comparisons were performed, and IBM SPSS Statistics 26 was used to process the data.

Comparative analysis of the chemical properties of soils of different vegetation types

The soil chemical properties of the four vegetation types were significantly different (p < 0.01; Table 1). pH value (mean 6.91), total soil carbon (mean 208.9 g/kg), total nitrogen(mean 13.69 g/kg), C/N (mean 15.27), N/P (mean 33.62) and C/P (mean 513.13) were all highest in DP and were significantly higher than in the other three vegetation types. The mean content of total phosphorus was 0.48 g/kg for ZM, which was significantly higher than the others (p = 2.61 × 10⁻³). HH and DP had significantly higher contents of each soil chemical property than CR and ZM.

Comparative analysis of soil microbial community diversity in different vegetation types

In total, a total of 1,143,966 high-quality bacterial sequences and 1,040,630 high-quality fungal sequences were obtained from all samples, which provided an opportunity to delve deeper into the bacterial and fungal communities. On average, each sample contained 71,498 bacterial sequences and 65,039 fungal sequences. The number of bacterial sequences ranged from 54,510 to 90,097 per sample, while the number of fungal sequences ranged from 51,968 to 77,522 per sample. All sequences were assigned to a total of 72,547 bacterial ASVs and 61,430 fungal ASVs. The number of soil bacterial ASVs in HH, DP, ZM, and CR was 10,481, 10,518, 10,226, 8,874, respectively (Fig. 1A). The number of ASVs in soil fungi was 1,410, 1,090, 563, 384 (Fig. 1B).

The number of shared soil bacterial ASVs among HH, DP, ZM, and CR was 119 (Fig. 1A). The 119 ASVs belonged to the phyla Acidobacteria, Actinobacteria, Chloroflexi, Firmicutes Gemmatimonadetes, Nitrospirae, Proteobacteria, Rokubacteria, Verrucomicrobia (Fig. 2A). Among them, Proteobacteria and Actinobacteria accounted for the highest percentage, with 45.39% and 33.61% respectively. Alphaproteobacteria was the dominant order shared by the four samples, with an ASV count of 46. The number of shared soil fungal ASVs among HH, DP, ZM, and CR was 11 (Fig. 1B). The ASVs common to all four samples were in the phyla Ascomycota (ASVs 81.82%) and Basidiomycota (ASVs 18.18%) (Fig. 2B). Dothideomycetes and Sordariomycetes were the dominant orders shared by the four samples.

The soil alpha diversity analysis of the different vegetation types was shown in Fig. 3. As a whole, the Chao 1 index, Shannon index and Pielou’s evenness index were significantly different between the bacterial communities of the four soil samples (Fig. 3A). HH had the highest Chao 1 index (mean 4572.95), Shannon index (mean 10.88) and Pielou’s evenness index (mean 0.91), while the lowest mean Chao 1 index, Shannon index and Pielou’s evenness index were found in CR at 3628.18, 10.22 and 0.88, respectively. Of the four soil samples, only HH and CR differed significantly between the two pairs (Chao 1 index p = 0.014; Shannon index p = 0.0065 and Pielou index p = 0.005; Table S1). The rest of the treatments were not significantly different between the two pairs (Table S1). Unlike the bacterial community, only the Chao 1 index and Shannon index were significantly different between the fungal communities of the four soil samples, while the Pielou’s evenness index was not significantly different between them (Fig. 3B). The mean values of the Chao 1 index and Shannon index followed the same pattern as the bacterial community, being HH (527.34; 6.87) >DP (489.56; 6.34) >ZM (231.67; 6.06) >CR (155.71; 5.54). The Chao 1 index and Shannon index of HH were significantly higher than that of CR (p = 0.0065; Table S1). The Chao 1 index of CR was significantly lower than that of DP (p = 0.038; Table S1). In contrast to the bacterial community, ZM had the highest mean Pielou’s evenness index and DP the lowest. the mean Pielou’s evenness index of HH was only 0.0004 higher than that of CR. Soil bacterial communities had higher Chao 1 index, Shannon index, and Pielou’s evenness index than soil fungal communities.

Comparative analysis of soil microbial community composition in different vegetation types

The relative abundance of soil microbial (bacterial and fungal) communities in each of the four vegetation types was analyzed at the phylum level and genus level, and communities with relative abundances greater than 1% were selected for comparative analysis (Table S2). Eight bacterial phyla with relative abundance greater than 1% were Actinobacteria, Proteobacteria, Acidobacteria, Chloroflexi, Gemmatimonadetes, Rokubacteria, Planctomycetes, Bacteroidetes (Fig. 4A). Among the four vegetation types, the phylum Actinobacteria, Proteobacteria, Acidobacteria, and Chloroflexi had average relative abundances greater than 10%, and they were also among the top 5 phyla in terms of the number of ASVs shared by all soil bacterial communities (Fig. 2A, Table S2). HH had the highest relative abundance of Proteobacteria and Acidobacteria of all samples at 36.31% and 15.38%, while Actinobacteria had the lowest relative abundance at 24.81%. CR had the highest relative abundance of Actinobacteria (42.99%) and Chloroflexi (18.57%) and the lowest of Proteobacteria (17.39%). Acidobacteria had the lowest relative abundance of 9.63% in ZM, while Chloroflexi had the lowest relative abundance of 3.28% in DP, the lowest of all.

At the fungal phylum level, relative abundances greater than 1% were Ascomycota, Basidiomycota, Mortierellomycota (Table S2). The only soil fungal phyla with average relative abundances greater than 10% were Ascomycota and Basidiomycota, which were also the phyla to which all the sampled fungal communities shared ASVs (Fig. 2B, Table S2). The relative abundance of the remaining eight phyla, with the exception of Mortierellomycota, was not higher than 0.01%. CR had the highest relative abundance of Ascomycota (85.46%) and the lowest relative abundance of Basidiomycota (12.44%). However, in DP, Ascomycota had the lowest relative abundance at 55.17% and Basidiomycota had the highest relative abundance at 31.91%.

At the level of soil bacterial genera in the four vegetation types, CR had not only the highest relative abundance of AD3 (7.62%), Acidothermus (7.70%), Rubrobacter (2.16%) and JG30-KF-AS9 (2.87%), but also the lowest relative abundance of Subgroup_6 (1.97%), 67-14 (1.42%), Rokubacteriales (1.40%), bacteriap25 (0.65%), Gaiella (1.61%), Solirubrobacter (0.87%). HH had the highest relative abundance of Subgroup_6 (8.92%), Rokubacteriales (4.44%), bacteriap25 (4.64%), Gaiella (1.90%) and the lowest relative abundance of AD3 (0.001%), Rubrobacter (0.82%) and JG30-KF-AS9 (0.01%) had the lowest relative abundance. The highest relative abundance of DP’s 67-14 (7.25%) and Solirubrobacter (2.30%) was found among them, while Acidothermus had the lowest relative abundance of 0.007% (Table S2). At the soil fungal genus level, the top ten genera in terms of CR relative abundance totaled 54.18%, the highest of them all, with ZM in second place at only 27.01% and HH the lowest at 21.11% (Table S2). Fusarium had the highest relative abundance in CR(17.84%) and the lowest relative abundance in DP (1.997%). Aspergillus had the highest relative abundance in HH (5.89%) and the lowest in DP (3.34%). Hygrocybe had a relative abundance of 17.03% in DP but was not distributed in HH and ZM. Paramyrothecium was not contained in DP and ZM and Acrophialophora was not distributed in HH and DP.

Based on the soil microbial community order level, the four sample sites were clustered using the average algorithm. At the bacterial order level, the HH and DP communities were combined into a single branch, and the ZM and CR communities were in the same branch (Fig. 4A). The bacterial communities of the two sample sites in the same branch, i.e., between HH and DP and between ZM and CR, were more similar. At the fungal order level, the ZM and CR communities merged into one branch and then coalesced with DP; the HH community was a separate branch (Fig. 4B). This indicated that the CR and ZM communities were more similar and least similar to HH. This clustering result reflected that the soil microbial community characteristics of the four sample site communities were closely related to the sample site plant community composition.

Further, a PERMANOVA test was conducted to examine the differences in soil microbial community composition between the four vegetation types (Table 2). As a result, the differences in soil microbial community composition between the four samples reached a significant level (p < 0.05). Among them, DP and ZM had the most significant differences in microbial community composition. The study showed that the type of above-ground vegetation was an important driver of structural differences in soil microbial community composition, and that differences in soil microbial community composition were more pronounced between vegetation types.

Relationship between soil nutrient factors and microbial community composition

The relative abundance of microphyla, with the exception of Basidiomycota, showed significant correlations with soil nutrient factors when their relative abundances exceeded 1% (Fig. 5). Specifically, the relative abundances of Actinobacteria, Chloroflexi, Gemmatimonadetes, and Ascomycota were significantly negatively correlated with soil pH (r =  − 0.55, p < 0.05; r =  − 0.68, p < 0.01; r =  − 0.72, p < 0.01; r =  − 0.52, p < 0.05), total carbon (r =  − 0.68, p < 0.01; r =  − 0.89, p < 0.01; r =  − 0.75, p < 0.01; r =  − 0.75, p < 0.01), total nitrogen (r =  − 0.68, p < 0.01; r =  − 0.89, p < 0.01; r =  − 0.74, p < 0.01; r =  − 0.71, p < 0.01), C/N (r =  − 0.65, p < 0.01; r =  − 0.77, p < 0.01; r =  − 0.94, p < 0.01; r =  − 0.62, p < 0.05), N/P (r =  − 0.66, p < 0.01; r =  − 0.77, p < 0.01; r =  − 0.94, p < 0.01; r =  − 0.62, p < 0.01), and C/P (r =  − 0.66, p < 0.01; r =  − 0.91, p < 0.01; r =  − 0.77, p < 0.01; r =  − 0.71, p < 0.01). On the other hand, the relative abundances of Acidobacteria, Rokubacteria, and Planctomycetes were significantly positively correlated with soil pH (r = 0.53, p < 0.05; r = 0.50, p < 0.05; r = 0.63, p < 0.01), total carbon (r = 0.64, p < 0.01; r = 0.73, p < 0.01; r = 0.71, p < 0.01), total nitrogen (r = 0.63, p < 0.01; r = 0.72, p < 0.01; r = 0.71, p < 0.01), C/N (r = 0.83, p < 0.01; r = 0.62, p < 0.05; r = 0.68, p < 0.01), N/P (r = 0.63, p < 0.01; r = 0.77, p < 0.01; r = 0.72, p < 0.01), and C/P (r = 0.64, p < 0.01; r = 0.76, p < 0.01; r = 0.70, p < 0.01). Furthermore, the relative abundances of Proteobacteria and Bacteroidetes were significantly positively correlated with soil total carbon (r = 0.73, p < 0.01; r = 0.78, p < 0.01), total nitrogen (r = 0.73, p < 0.01; r = 0.78, p < 0.01), C/N (r = 0.56, p < 0.05; r = 0.57, p < 0.05), N/P (r = 0.73, p < 0.01; r = 0.79, p < 0.01), and C/P (r = 0.73, p < 0.01; r = 0.79, p < 0.01), while showing significant negative correlations with soil total phosphorus (r =  − 0.52, p < 0.05; r = −0.61, p < 0.05). Additionally, the relative abundance of Rokubacteria and Mortierellomycota were negatively correlated with soil total phosphorus content (r =  − 0.55, p < 0.05; r =  − 0.57, p < 0.05). Notably, there was no significant correlation observed between the relative abundance of Basidiomycota and soil nutrient content (p > 0.05). Among these microphyla, Gemmatimonadetes showed the strongest negative correlation with soil pH (r =  − 0.72) and C/N (r =  − 0.94), Chloroflexi had the highest negative correlation with soil total carbon (r =  − 0.89), total nitrogen (r =  − 0.89), N/P (r =  − 0.92), and C/P (r =  − 0.91), and Bacteroidetes had the highest negative correlation with soil total phosphorus content (r =  − 0.61).

At the genus level, our analysis revealed significant associations between soil nutrient factors and bacterial genera, with relative abundance exceeding 1% (Fig. 6A). Soil pH exhibited significant correlations with the relative abundance of AD3, Acidothermus, Rokubacteriales, Gaiella, JG30-KF-AS9, and IMCC26256 (p < 0.05). Among these genera, only Rokubacteriales showed a positive correlation (r > 0). Notably, the relative abundances of Subgroup_6, 67-14, Rokubacteriales, bacteriap25, and Solirubrobacter displayed significant positive correlations with total carbon, total nitrogen, N/P, and C/P ratios (r > 0, p < 0.05). Conversely, these nutrient factors exhibited negative correlations with the abundance of AD3, Acidothermus, JG30-KF-AS9, IMCC26256, and TK10 (r <0, p < 0.05). Additionally, soil C/N ratio demonstrated significant positive correlations with the relative abundance of Subgroup_6, Rokubacteriales, and bacteriap25 (r > 0, p < 0.05), while exhibiting negative correlations with the abundance of AD3, Acidothermus, JG30-KF-AS9, IMCC26256, and TK10 (r < 0, p < 0.05). Furthermore, only Subgroup_6, Rokubacteriales, and Gaiella displayed negative correlations with soil total phosphorus content (r < 0, p < 0.05), whereas the relative abundances of Rubrobacter, TK10, and RB41 exhibited positive correlations with soil total phosphorus content (r > 0, p < 0.05). Notably, the relative abundances of Mycobacterium and Bradyrhizobium did not exhibit any significant correlations with the examined soil nutrient factors (p > 0.05). The relative abundances of Fusarium, Talaromyces, Basidioscus, Acrophialophora and Dokmaia exhibited negative correlations with the other factors, with the exception of soil total phosphorus (r < 0, p < 0.05). Conversely, a significant positive correlation was observed with the relative abundance of Hygrocybe (r > 0, p < 0.05) (Fig. 6B). The relative abundances of Penicillium, Humicola and Saitozyma were negatively correlated with soil total carbon, total nitrogen, N/P and C/P (r < 0, p < 0.05). Additionally, the relative abundance of Saitozyma displayed a negative correlation with soil C/N (r < 0, p < 0.01). The relative abundances of Aspergillus and Mortierella were negatively correlated solely with soil total phosphorus content (r < 0, p < 0.05). On the other hand, the relative abundance of Geastrum exhibited significant positive correlations with soil total carbon, total nitrogen, C/N, N/P and C/P (r > 0, p < 0.05), while displaying a negative correlation with soil total phosphorus content (r < 0, p < 0.05).