MATERIALS AND METHODS
2.1 Site selection and soil collection
Field sampling was conducted in Guizhou Province and the Guangxi Autonomous Region of southwestern China. The spatial range was 104°82′–108°37′E, 22°42′–27°53′N (Fig. 1). The soil type was calcareous lithosols (limestone soil) according to the Soil Taxonomy of China and the FAO/UNESCO system of soil classification. Samples were collected across six counties from typical karst regions with plantation forest (PF), shrubland (SH), and reference cropland (CR). Plantation forest and shrubland comprised artificial and natural vegetation restoration, respectively, and cultivation was abandoned there in 2002–2003. Shuicheng, Jinsha, and Duyun counties are in Guizhou Province and Huanjiang, Mashan, and Longzhou counties are in Guangxi Autonomous Region. The mean annual temperatures (MAT) within the Guizhou and Guangxi sampling sites were 14.6 °C and 20.9 °C, respectively. The mean annual precipitation (MAP) levels in Guizhou and Guangxi were 1,100 mm and 1,260 mm, respectively. Further details of the climate conditions, geographic locations, and dominant plant species are listed in Table S1.
Soil sample collection was performed between August and October 2018. There were 54 soil samples (three land use types × six sampling counties × three repetitions). In each 30 m × 30 m sampling field, twenty 38-mm-diameter soil cores were collected at 15 cm soil depth and pooled. Each soil sample was passed through a 2-mm sieve to remove rocks, roots, and organic debris. Physicochemical properties were determined for air-dried soil. Genomic DNA was extracted from soil stored at -80 °C.
2 .2 Soil physicochemical property measurements
The pH of 10 g soil suspended in 25 mL water was measured with a glass electrode pH meter. Exchangeable calcium (Ca) and magnesium (Mg) were displaced via compulsive exchange in 1 M ammonium acetate (pH 7). Soil organic carbon (SOC) was measured by dichromate redox colorimetry. Particulate organic C (POC) and mineral-associated organic C (MOC) were isolated by size fractionation (> 53 μm and < 53 μm) and analyzed in the same way as the bulk SOC. Total N (TN) was determined using an elemental (CN) analyzer (Vario MAX CN; Elementar, Hanau, Germany). Ammonium N (NH4+) and nitrate N (NO3) were extracted from 10 g fresh soil by extraction with 0.5 M KCl and measured with an auto flow analyzer. Total phosphorus (TP) and available phosphorus (AP) were determined by digestion in a solution containing H2SO4 and HClO4 and by molybdenum blue colorimetry, respectively (Hu et al., 2021; Xiao et al., 2019, 2020).
2.3 Soil DNAextraction and amplicon sequencing
The total DNA was extracted from 0.5 g soil samples using a FastDNA Spin Kit (MP Biomedicals, Santa Ana, CA, USA) according to the manufacturer’s instructions. DNA quality and quantity were determined with a NanoDrop 2000 UV-vis spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). Bacterial 16S rDNA gene was amplified using the primer pairs 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’). Fungal ITS and protist 18S rRNA genes were amplified with the primer sets ITS1F (5’-CTTGGTCATTTAGAGGAAGTAA-3’)/ITS2R (5’-GCTGCGTTCTTCATCGATGC-3’) and TAReuk454FWD1 (5’-CCAGCA(G/C)C(C/T)GCGGTAATTCC-3’)/TAReukREV3 (5’-ACTTTCGTTCTTGAT (C/T)(A/G)A-3’), respectively. Samples were distinguished by attaching barcode sequences to each primer.
PCR was performed in a 20 µL mixture containing 4 µL of 5× FastPfu Buffer, 2 µL of 2.5 mM dNTPs, 0.8 µL of 5 μM forward and reverse primers, 0.4 µL TaKaRa rTaq DNA polymerase (TaKaRa Bio Inc., Shiga, Japan), 0.2 µL BSA, 10 ng template DNA, and 20 µL ddH2O. The thermal cycles for the 16S rDNA and ITS genes were, 95 °C for 3 min, 27 cycles and 35 cycles for 16S rDNA and ITS genes, respectively, 95 °C for 30 s, 55 °C for 30 s, 72 °C for 45 s, and a final extension at 72 °C for 10 min. The PCR program used for protist 18S rRNA gene was, initial denaturation at 95 °C for 5 min, 10 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C for 45 s, extension at 72 °C for 60 s, 25 cycles of denaturation at 94 °C for 30 s, annealing at 45 °C, 47 °C, 48 °C, and 49 °C for 45 s each, extension at 72 °C for 60 s, and a final extension at 72 °C for 2 min. PCR was performed at Shanghai Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China) on an Illumina MiSeq PE300 platform (Illumina, San Diego, CA, USA).
2.4 Sequence analysis
Bioinformatic analyses of all three amplicons were performed with QIIME 2 2021.2 (Bolyen et al., 2019). The raw reads were imported into the QIIME 2 environment and denoised with the DADA2 algorithm via the q2-data2 plugin (Callahan et al., 2016). Operational taxonomic units (OTUs) were selected at 97% sequence similarity in the q2-vsearch plugin (Rognes et al., 2016). OTU taxonomies for bacteria, fungi, and protist were assigned using the q2-feature-classifier (Bokulich et al., 2018) vs. the SILVA_v132 (Quast et al., 2012), UNITE_v8.3 (Nilsson et al., 2019), and PR2 (V4.13.0) (Guillou et al., 2012) reference sequences, respectively. Alpha-diversity metrics (richness and Shannon index) were calculated in the “RAM” package of R (R Core Team, Vienna, Austria).
2.5 Statistical analysis
Data were processed in R v. 3.6.3 (R Core Team). Significant differences in microbial diversity, microbial groups, and environmental factors among vegetation types were analyzed by one-way ANOVA and Duncan’s post hoc test (p < 0.05). Three-way nested PERMANOVA was used to determine the main and interactive effects of MAT, MAP, and vegetation type on microbial diversity and community composition (Xiao et al., 2019). Non-metric multidimensional scaling (NMDS) was performed to reveal dissimilarities in microbial, fungal, and protist community composition among treatments. Pearson’s correlation coefficients and the Mantel test were applied to determine the relationships among environmental factors (geographical factors, climatic conditions, and soil properties), microbial diversity, and community composition.
To determine the structure and identify interactions among microbial groups, bacterial, fungal, and protist community networks were constructed across all samples and vegetation types under low (GZ) and high (GX) climate conditions. Spearman correlations were calculated to clarify the associations among the microbial taxa detected in > 1/3 of all samples per network. Significant Spearman associations (ρ > 0.75; 1,000-fold permutations; p< 0.05) and the corresponding OTUs served as network nodes and edges, respectively. Positive and negative associations suggested cooperative and competitive relationships, respectively. Network properties such as node, edge, density, and degree were calculated in each network using the “igraph” package in R (R Core Team) (Csardi and Nepusz, 2006).
Structural equation modeling (SEM) based on SPSS Amos 21 (IBM Corp., Armonk, NY, USA) was conducted to evaluate the causal relationships among MAT, vegetation type, soil property, and microbial community composition (Hu et al., 2021; Yang et al., 2017). To simplify the SEM, a principal component analysis (PCA) of the first component (PC1) was applied to represent soil properties and microbial community compositions based on bacterial, fungal, and protist NMDS axes 1 and 2. Soil properties included pH, soil exchangeable Ca2+ (Ca), soil exchangeable Mg2+(Mg), soil organic carbon (SOC), particulate organic carbon (POC), mineral-associated organic carbon (MOC), total nitrogen (TN), ratio of soil carbon to soil nitrogen (C/N), ammonium nitrogen (NH4+), nitrate nitrogen (NO3), total phosphorus (TP), and available phosphorus (AP).
3 RESULTS
3.1Bacterial, fungal, and protist diversity and community composition
Only bacterial diversity was affected by temperature (Table 1). Bacterial richness and the Shannon indices for plantation forest and shrubland were lower in Guangxi (high temperature) than Guizhou (low temperature) (Fig. 2a). Vegetation type significantly affected bacterial, fungal, and protist diversity (Table 1). Bacterial richness and Shannon index were lower under shrubland than cropland (except for richness in Guizhou) (Fig. 2a). Fungal richness was higher in shrubland than cropland under both climate conditions (Fig. 2b). At Guangxi, protist Shannon index was lower under plantation forest and shrubland than cropland. Cropland protist richness was higher than that of shrubland at Guangxi (Fig. 2c). There were no differences between plantation forest and shrubland in terms of bacterial, fungal, and or protistan diversity (Fig. 2).
Across all sites, phylum-level bacterial community composition was dominated by Actinobacteria (29.3%), Proteobacteria (25.6%), and Acidobacteria (16.7%) (Fig. 3a). The most abundant fungal phyla were Ascomycota (45.0%) and Basidiomycota (27.9%) (Fig. 3b). The main protist phyla were Cercozoa (33.9%) and Apicomplexa (22.1%) (Fig. 3c). NMDS ordination revealed that bacterial, fungal, and protist community composition differed among vegetation types between Guizhou and Guangxi. Bacterial, fungal, and protist community composition differed between cropland and vegetation recovery but were similar between plantation forest and shrubland (Figs. 3d-f). Temperature and vegetation type significantly affected bacterial, fungal, and protist community composition (Table 1).
The relative abundance of Proteobacteria (bacteria) under all vegetation types was lower at Guangxi than Guizhou (Table S2, Fig. S2). The relative abundance of Mucoromycota (fungi) was higher under plantation forest and shrubland at Guizhou than it was at Guangxi (Table S3, Fig. S3). The relative abundance of Cercozoa (protist) was higher under plantation forest and shrubland at Guizhou than it was at Guangxi (Table S4, Fig. S4). Relative bacterial phyla Cyanobacteria, Gemmatimonadetes, and Nitrospirae abundances were higher under cropland than plantation forest or shrubland. Relative Proteobacteria and Verrucomicrobia abundances were higher under shrubland than cropland (Table S2, Fig. S2). Relative fungal phyla Ascomycota and Mucoromycota abundances were higher under cropland than plantation forest and shrubland. Nevertheless, the opposite was true for Basidiomycota (Table S3, Fig. S3). Relative protist phyla Ciliophora, Lobosa, and Ochrophyta abundances were higher under cropland than plantation forest or shrubland. However, the opposite was true for Apicomplexa (Table S4, Fig. S4).
3.2 Co-occurrence networks amongbacterial, fungal, and protist taxa
Microbiome network complexity was greater at Guizhou (GZ) than Guangxi (GX). There were more nodes and edges under GZ than GX. The most densely connected network occurred under cropland (Table 2, Fig. 4). There were fewer positive and negative correlation links among bacterial, fungal, and protist taxa under vegetation restoration than cropland. Shrubland had the fewest positive and negative correlation links among microbial taxa (Fig. 4). Node connectivity was greater under cropland than plantation forest or shrubland. Shrubland had the fewest node connections. Under all vegetation and both climate types, bacteria had the most nodes, followed by fungi and protists (Table 2, Fig. S5).
Acidobacteria (bacteria), Ascomycota (fungi), and Cercozoa (protist) formed remarkably centralized networks and showed the most links with other microbial groups (Table S5, Fig. 4). Bacterial phyla Acidobacteria and Proteobacteria were the most closely connected with Ascomycota and Cercozoa. The strongest and most numerous links between fungal and protist taxa were determined for Ascomycota and Cercozoa, followed by Ascomycota and Lobosa (Table S5). There were more interactions between bacterial and fungal taxa (edge links) than there were between bacterial and protist taxa or between fungal and protist taxa. The proportions of edges linking bacteria to protists were lower under plantation forest and shrubland than cropland. Conversely, the proportions of edges linking fungi to protists were higher under plantation forest and shrubland than cropland (Table 2).
3.3 Relationships among environmental factors and bacterial, fungal, and protist communities
The various bacterial, fungal, and protist taxa had different effects on soil carbon (SOC, POC, and MOC), nitrogen (TN, NH4+, and NO3), and phosphorus (TP and AP). More bacterial, fungal, and protist OTUs correlated with C and N than they did with P (Fig. 5). For bacterial taxa, soil SOC, POC, and MOC were most strongly correlated with Proteobacteria. Acidobacteria and Proteobacteria had the most OTUs correlated with soil TN, NH4+, NO3, TP, and AP (Fig. 5a). Soil carbon (SOC, POC, and MOC), nitrogen (TN, NH4+, and NO3), and phosphorus (TP and AP) were strongly correlated with fungal phyla Ascomycota (Fig. 5b). Protist phyla Cercozoa were strongly correlated with soil C, N, and P (Fig. 5c).
Pearson’s correlation coefficients showed that bacterial richness and Shannon index were negatively correlated with MAT, MAP, and MOC but positively correlated with latitude and altitude (Fig. 6a). Fungal diversity was not significantly correlated with any geographical factor, climate condition, or soil property (Fig. 6b). Protist diversity was positively correlated with AP but negatively correlated with SOC, POC, TN, and NH4+ (Fig. 6c). Bacterial, fungal, and protist community compositions were strongly affected by environmental factors (geography, climate, and soil physicochemistry) (Fig. 6).
Structural equation modeling (SEM) evaluated the linkages among environmental factors and microbial community composition, and ranked the relative importance of each driver on the microbial communities. Temperature and vegetation type directly influenced soil property and indirectly affected bacterial, fungal, and protist community composition. MAT had greater standard total effect on bacterial and fungal community composition than it did on protist community composition. Vegetation type had the highest standard total effect on protist community composition followed by fungal community composition and bacterial community composition. Bacterial and fungal community composition was affected by MAT and vegetation type. Moreover, vegetation type and soil property had greater influences than MAT on protist community composition (Fig. 7).