Soil pH Is the Primary Factor Correlating With Soil Microbiome in Karst Rocky Desertification Regions in the Wushan County, Chongqing, China.

Karst rocky desertification (KRD) is a process of land degradation, which causes desert-like landscapes, deconstruction of endemic biomass, and declined soil quality. The relationship of KRD progression with above-ground communities (e.g. vegetation and animal) is well-studied. Interaction of soil desertification with underground communities, such as soil microbiome, however, is vastly unknown. This study characterizes change in soil bacterial community in response to KRD progression. Soil bacterial communities were surveyed by deep sequencing of 16S amplicons. Eight soil properties, pH, soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK), were measured to assess soil quality. We find that the overall soil quality decreases along with KRD progressive gradient. Soil bacterial community compositions are distinguishingly different in KRD stages. The richness and diversity in bacterial community do not significantly change with KRD progression although a slight increase in diversity was observed. A slight decrease in richness was seen in SKRD areas. Soil pH primarily correlates with bacterial community composition. We identified a core microbiome for KRD soils consisting of; Acidobacteria, Alpha-Proteobacteria, Planctomycetes, Beta-Proteobacteria, Actinobacteria, Firmicutes, Delta-Proteobacteria, Chloroflexi, Bacteroidetes, Nitrospirae, and Gemmatimonadetes in this study. Phylum Cyanobacteria is significantly abundant in non-degraded soils, suggesting that Cyanobacterial activities might be correlated to soil quality. Our results suggest that Proteobacteria are sensitive to changes in soil properties caused by the KRD progression. Alpha- and beta-Proteobacteria significantly predominated in SKRD compared to NKRD, suggesting that Proteobacteria, along with many others in the core microbiome (Acidobacteria, Actinobacteria, Firmicutes, and Nitrospirae), were active in nutrient limiting degraded soils. This study demonstrates the relationship of soil properties with bacterial community in KRD areas. Our results fill the gap of knowledge on change in soil bacterial community during KRD progression.

Introduction

Karst is a type of natural landscape with a fragile ecosystem. Land in karst area is highly susceptible to envn class="Chemical">ironmental disturbances (Calò and Parise, 2006). Soil structure in the karst region is unstable and prone to deterioration (Gutiérrez et al., 2014). Karst rocky desertification (KRD) is the progress of soil desertification. This kind of desertification often occurs in the karst area, resulting in desert-like landscape (Yan and Cai, 2015). Increased soil erosion and bedrock exposure along with decreased soil quality and biomass are main consequences of KRD progression (Ford and Williams, 1989). Decreased soil quality is defined as low biological productivity, environmental quality, and plant/animal health (Karlen et al., 1997). Severe soil desertification deprives cultivated land and affects global geomorphology (Jiang et al., 2014). Altogether, the progression of KRD has been drawing considerable attention. Karst landscape widely exists in southwest of China, especially in Yunan, Guizhou, and Chongqing (Peng et al., 2011). Karst area is classified into no-, light-, moderate-, and severe-KRD (NKRD, LKRD, MKRD, and SKRD); based on bedrock exposure rate, vegetation, and soil cover rate (Li et al., 2009). The land in the Wushan County (northeast Chongqing) is highly susceptible to desertification. Roughly 70% of soils (2066 km2) in the Wushan County are M- and S-KRD areas (Qi et al., 2017).

Many studies report causes of KRD and changes in soil and vegetation community in response to KRD progression (Calò and Parise, 2006; Xie et al., 2015; Zhang et al., 2016; Qi et al., 2017). It is suggested that besides the natural susceptibility to desertification in the karst region, anthropogenic disturbances and soil exploitation are the other main reasons accelerating soil desertification (Yan and Cai, 2015). Furthermore, soil pH significantly increases along the KRD gradient, and is associated with decreased soil quality (Wang et al., 2004). As the result, vegetation diversity and coverage significantly decreases in accordance to desertification progression (Zeng et al., 2007; Qi et al., 2017). Many studies focus on determining correlations of aboveground factors with KRD progression. Few studies, however, aim to understand change in underground community, such as soil microbiome, influenced by the KRD progression.

Soil microbiome is highly heterogeneous. Microbiome includes all the microorganisms that belong to bacteria, fungi, and archaea. n class="Species">Soil microbiome, therefore, is complex in structure and enormous in number. More importantly, soil microbiome has large impact on plant health, soil productivity, and ecosystem (Chaparro et al., 2012; Zolla et al., 2013; Lakshmanan et al., 2014). Soil bacteria play an essential role in soil weathering and nutrient cycling (Uroz et al., 2011). For instance, potassium solubilizing bacteria (KSB) produce low molecular weight organic acids that result in rock weathering; while other bacteria such as ammonia oxidizing bacteria (AOB) and nitrite oxidizing bacteria (NOB), produce enzymes that oxidize nitrogen compounds, thereby contributing to nutrient cycling (De Boer and Kowalchuk, 2001; Kowalchuk and Stephen, 2001; Mursyida et al., 2015). Copiotrophic bacteria, such as Firmicutes, Proteobacteria, and Actinobacteria, usually claim the predominance in fertile soils (Yao et al., 2000; Fierer et al., 2007). Oligotrophic taxa, such as Acidobacteria, on the other hand, are abundant in barren soils (Shah et al., 2011; Rime et al., 2015). Soil bacteria are particularly essential in terrestrial ecosystem, because (1) terrestrial soil is usually not supplemented with external fertilizer and (2) soil bacterial activities and byproducts affect soil fertility in the terrestrial ecosystem. It is important to characterize change in soil microbiome in response to KRD progression. Establishing changes in microbiome as well as the core microbiome helps to understand underground community composition and microbial contributions on the ecosystem in KRD areas.

In this study, we characterize change in soil bacterial community in accordance to KRD progression. Soil bacterial community was surveyed by deep sequencing of 16S hypervariable 4-5 region (V4-5) and Illumina HiSeq platform. Eight soil properties, pH, soil organic matter (SOM), total and available n class="Chemical">nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK), were measured for soil quality. We hypothesized that (1) the composition of soil bacterial community changed along the KRD gradient; (2) changes in soil properties caused by desertification impacted on soil bacterial community; and (3) certain bacteria were sensitive to KRD progression. To address our hypotheses, we (1) surveyed soil bacterial community and identified core microbiome in KRD regions; (2) determined soil property as the primary influential factor shaping bacterial community; and (3) analyzed correlation of soil properties with core phyla found in this study.

Materials and methods

Study area

Detailed descriptions of study area are reported in a previous study (Qi et al., 2017). Briefly, the study area is located in the Wushan County, northeast Chongqing (N 30°46′-31°28′, E 109°33′-110°11′). The main soil types in study area are Xanthic Ferralsols and Haplic n class="Chemical">Luvisols. Soil types are classified according to guidelines of Food and Agriculture Organization of the United Nations (FAO) (Gong, 2001). Xanthic Ferralsols is generally found in the hilly and low-mid mountain zone that is below 1,500 m of elevation. The Haplic Luvisols is mainly found in degraded soils with anticlinal structure and is neutral to slightly alkaline soil pH. The bedrock in study site is composed of limestone, clastic carbonate, and dolomite.

Study sites were evaluated and classified into four groups, no KRD (NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD), based on vegetation coverage rate and bedrock exposure rate (Table 1). A total of 12 sampling quadrats (20 × 20 m) were chosen, three for each stage of KRD. Three sub-sampling quadrats (10 × 10 m) were uniformly chosen within each sampling quadrat in order to study n class="Species">soil microbiome community comprehensively.

Table 1

Features and evaluations of karst rocky desertification areas.

Study sites	NKRD	LKRD	MKRD	SKRD
Coordinates	N 31°8′4″ E 110°0′45″N 31°9′1″ E 110°1′31″N 31°8′48″ E 110°1′40″N 31°7′54″ E 110°0′6″N 31°2′7″ E 110°44′11″	N 31°3′32″ E 109°45′56″N 31°4′45″ E 109°54′28″N 31°0′4″ E 109°51′8″N 31°0′3″ E 109°51′22″N 31°4′21″ E 109°48′59″	N 31°50′56″ E 110°2′41″N 31°0′1″ E 109°50′56″N 31°0′26″ E 109°52′16″N 31°9′30″ E 109°46′14″N 31°7′27″ E 109°45′24″	N 31°4′42″ E 109°54′25″N 31°2′20″ E 109°44′56″N 31°0′15″ E 109°51′53″N 31°1′55″ E 109°54′6″N 31°6′49″ E 109°44′34″
Elevation (M ASL)	850–1,100	750–1,100	850–1,200	900–1,200
Vegetation	Coniferous forest	Shrub-land	Shrub-grassland	Shrub-grassland
Vegetation Coverage (%)	≥80	≥80	40–60	15–30
Bedrock	Limestone	Limestone	Limestone	Limestone
Bedrock Exposure Rate (%)	0	0–30	50–70	>70
Soil Type	Xanthic Ferralsols	Haplic Luvisols	Haplic Luvisols	Haplic Luvisols
Soil Depth (Cm)	40–70	20–40	20–30	10–20

Geographic features, soil type, bedrock exposure rate, and vegetation coverage are summarized in the table. Classification of KRD is based on vegetation coverage, soil depth, and bedrock exposure rates. No KRD (NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD).

Geographic features, soil type, bedrock exposure rate, and vegetation coverage are summarized in the table. Classification of KRD is based on vegetation coverage, soil depth, and bedrock exposure rates. No KRD (n class="Chemical">NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD).

Soil sample collection

Soil samples were collected in N-, L-, M-, and S-KRD regions with little n class="Species">human interference (i.e., land with no agricultural activities, vegetation restoration, or soil exploitation etc.). All the collecting tools were autoclaved. Three soil samples (5 g) were randomly and uniformly collected from the top soil layer (0–10 cm) in each sub-sampling quadrat. Any plant residues, such as roots and leaves, were either carefully avoided or removed during soil sample collection. Soil samples were collected and stored in 50 mL sterile tubes on ice (approximately 0°C). Three samples from one sub-sampling quadrat were pooled into the same 50 mL tube (15 g) to comprehensively represent soil microbiome community in a quadrat. All pooled samples from sampling quadrats were labeled as NKRD_1, NKRD_2, NKRD_3, LKRD_1, LKRD_2, LKRD_3, MKRD_1, MKRD_2, MKRD_3, SKRD_1, SKRD_2, and SKRD_3. Once arrived in the lab, soil samples were immediately mixed by Large Capacity Mixer (Glas-Col, Terre Haute, IN, USA) at 1,600 rpm, 4°C for 30 min. Two sets of soil subsamples (2 g) were collected in a sterile environment. One set of samples was sent to the Key Laboratory of Eco-environments in Three Gorges Reservoir Region for the soil properties measurement. The other set of samples was used for DNA isolation.

Soil properties measurement

Eight soil properties were measured in this study, namely pH, soil organic matter (SOM), total and available n class="Chemical">nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Detailed protocols were reported in previous studies (Li and Shao, 2006; Qi et al., 2017). Briefly, soil pH was measured using a pH meter (FE20, Mettler-Toledo Instruments, China) with a soil-to-water ratio of 1:2.5 (w/v). Soil organic carbon (SOC) was measured using potassium dichromate oxidation (Nelson and Sommers, 1982). Soil organic matter was calculated by multiplying SOC values by 1.724 (Walkley and Black, 1934). TN, AN, TK, AK, TP, and AP were measured using LAQUAtwin ion meter kits (Horiba Instruments, Singapore).

DNA isolation

DNA was isolated by using the MoBio PowerSoil DNA extraction kit and following manufacturer's instructions (MoBio Laboratories, Carlsbad, CA, USA). A total of 0.5 g sub-sampled soil was used. In order to prevent DNA shearing, an additional incubation step was added for 10 min at 65°C, followed by bead beating for 2 min (Lauber et al., 2009). DNA samples were investigated by using a 0.8% (wt/vol) low melting point agarose gel. High quality DNA was isolated and purified by using an agarose gel DNA purification kit (TaKaRa Bio USA, Inc., Mountain View, CA, USA). Purified DNA samples were quantified by using NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Samples were stored at −20°C till use.

PCR amplification and 16S V4-5 region sequencing

An aliquot of 50 ng DNA template was used for PCR amplification. The 16S V4-5 region was amplified using primer set: 515F (5′-GTGCCAGCMGCCGCGG-3′) and 907R (5′-CCGTCAATTCMTTTRAGTTT-3′) with barcodes. PCR amplification was performed by using Phusion® High-Fidelity PCR Master Mix (n class="Chemical">New England Biolabs, Ipswich, MA, USA). Quantification and qualification of PCR products were performed by mixing the same amount of 1X loading buffer that contained SYBR green with PCR products. Samples were loaded on 2% electrophoresis agarose gel for detection. Gel purification was performed using Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany).

TruSeq DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, CA, USA) was used to conduct sequencing library by following manufacturer's recommendations. Index codes were then added. The library quality was verified using the Qubit 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system. The library was sequenced on an Illumina HiSeq2500 instrument (Illumina, Inc., San Diego, CA, USA), 250 bp paired-end raw reads were generated.

Bioinformatic analysis

Mothur (v. 1.39.5) was used for data mining following mothur MiSeq standard operating procedure (Schloss et al., 2009). Briefly, “make.contigs” command was used to assemble paired-end reads. Barcode and primer sequences were removed. Command “screen.seqs” was used to trim the data with the criteria “maxambig = 0, maxhomop = 8, and maxlength = 380.” Contigs that were longer than 380 bp, or contained longer than 8 homopolymers, or contained undermined base were excluded. “trim.seqs,” “pre.cluster,” and “chimera.uchime,” were used to further process sequences. Uchime was used to remove chimeric sequences (Edgar et al., 2011). Sequences were then classified using the Bayesian classifier with the “classif.seqs” command. Silva non-redundant v123 database was used as the reference (Quast et al., 2013). Archaea, chloroplasts, and mitochondria were excluded by using command “remove.lineage” with criteria “taxon = Chloroplast-Mitochondria-unknown-Archaea-Eukaryota.” Command “cluster.split” was used to assign qualified sequences into operational taxonomic units (OTUs) based on at least 97% similarity (OTU0.03). Command “classify.otu” was used to classify all the sequences of OTU0.03 into taxonomic groups at the bootstrap threshold of 80%. Random sub-sampling with the least amount of sequences among all the samples was conducted to avoid sequencing bias. The similarities among communities based on communities' memberships were measured by using unifrac-based metrics generated with command “unifrac.weighted” (Lozupone et al., 2011). Raw sequencing data was submitted to Sequence Read Archive (SRA) in n class="Chemical">NCBI, access number SRP124255.

Statistical analysis

One-way analysis of variance (One-way ANOVA) and Pearson correlation analysis were performed by using R 3.3.2 statistical software (Team, 2013). Tukey's (honest significant difference) n class="Disease">HSD test was used to determine differences of soil properties caused by soil degradation. The change in soil properties were considered significantly different if P < 0.05. Correlations of soil properties with microbiome were determined by using Pearson correlation analysis. Soil properties and microbiome were considered significantly correlated if the P < 0.05.

Linear discriminant analysis (LDA) Effect Size (LEfSe) was performed to determine the soil microbiome that significantly responded to KRD progression (Segata et al., 2011). More specifically, the non-parametric factorial Kruskal-Wallis (KW) sum-rank test was first used to detect microbiome with significant abundant differences between KRD stages (P < 0.05). The unpaired Wilcoxon rank-sum test was used to compare the significant abundant differences among taxa under influences of soil desertification (P < 0.05). Linear Discriminant Analysis was applied to calculate the effective size of abundant differences. The LAD scores were normalized by log10.

Distance-based redundancy analysis (db-RDA) was conducted by using vegan package in R (Oksanen et al., 2014). The db-RDA was used to perform the direct gradient analysis. Constrains and accumulates explanatory variables are scattered on the ordination exes (RDA1 and RDA2). The two ordinations axes (RDA1 and RDA2) were constructed by constraining an entire set of explanatory variables, eight soil properties in this study. Mantel test at 999 permutations was used to determine impact of soil properties on bacterial community composition. Impact of soil properties on microbiome community was considered significant if P < 0.05. Weighted unifrac-distance was used to determine soil bacterial community structural similarity. Communities with similar compositions closely clustered, and vice versa.

The multivariate regression tree (MRT) was conducted by using the mvpart package in R (De'Ath, 2002). The MRT tree was used to reveal the relationship of relative abundance of core phyla with changes in soil properties caused by KRD progression.

Results

Change in soil properties along with desertification progression

Changing trends of soil properties in KRD regions are shown in Figures 1A,B. Detailed values of soil properties are shown in Supplementary Table 1. One-way ANOVA (Supplementary Table 2) followed by Tukey's n class="Disease">HSD test (Table 2) were used to determine differences of soil properties among KRD areas. Overall, soil pH, SOM, and AK significantly increased along with degradation gradients (P < 0.05); whereas TK significantly decreased (P = 0.03). Total nitrogen in M- and S-KRD areas was significantly higher compared to that of in N- and L-KRD regions; whereas TP in SKRD was significantly lower than that of in NKRD regions.

Figure 1

Change in soil properties along with karst rocky desertification gradient (A,B). No KRD (NKRD) is in red color, latent KRD (LKRD) is in green color, moderate KRD (MKRD) is in yellow color, and severe KRD (SKRD) is in purple color. Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK).

Table 2

Average value of soil properties in karst rocky desertification regions.

	pH	SOM g/kg	TN g/kg	TP g/kg	TK g/kg	AN mg/g	AP mg/g	AK mg/g
NKRD	6.07c	10.36b	0.40b	0.32a	14.28ab	45.36a	3.03a	136.46b
LKRD	7.71b	27.63b	0.43ab	0.30ab	14.75a	47.77a	2.10a	146.44b
MKRD	8.05b	39.22ab	0.61a	0.24ab	10.72b	62.61a	2.44a	182.17a
SKRD	8.53a	47.47a	0.67a	0.18b	9.99b	65.97a	2.29a	178.89a

The Tukey's HSD test is used to determine the significant difference between variables. Different letters are assigned to significant difference, P < 0.05. No KRD (NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD). Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK).

Change in soil properties along with karst rocky desertification gradient (A,B). No KRD (n class="Chemical">NKRD) is in red color, latent KRD (LKRD) is in green color, moderate KRD (MKRD) is in yellow color, and severe KRD (SKRD) is in purple color. Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK).

The Tukey's HSD test is used to determine the significant difference between variables. Different letters are assigned to significant difference, P < 0.05. n class="Chemical">No KRD (NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD). Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK).

Soil microbiome in karst rocky desertification regions

Total of 12 soil samples were sequenced by targeting 16S V4-5 region, with three replicates in each KRD area. Soil samples were surveyed to study change in soil bacterial community along with karst desertification progression. A total of 1,340,457 paired-end raw reads were generated using Illumina HiSeq platform. Averages of 111,704 ± 28,034 paired-end raw reads were obtained per sample. A total of 1,181,535 qualified reads remained after data trimming and sub-sampling, with an average of 98,461 ± 25,371 reads per sample. A total of 25 phyla were classified with qualified reads. All the reads were further categorized at class, order, family, and genus levels. Overall, 54 classes, 116 orders, 205 families, and 487 genera were classified among all the samples. Detailed OTU information is shown in Supplementary Table 3.

Twelve phyla took up more than 1% of the total abundance, namely Acidobacteria (24.17%), Alpha-Proteobacteria (20.57%), Planctomycetes (11.51%), Beta-Proteobacteria (9.39%), Actinobacteria (6.90%), Firmicutes (5.97%), Deltaproteobacteria (4.80%), Chloroflexi (4.73%), Bacteroidetes (4.43%), Nitrospirae (2.13%), Gemmatimonadetes (1.29%), and Deinococcus-Thermus (1.26%). Core microbiome at phylum level was defined as phyla that were found in all 15 samples and took up at least 1% of the relative abundance within a sample. Under such criteria, phyla Acidobacteria, Alpha-Proteobacteria, Planctomycetes, Beta-Proteobacteria, Actinobacteria, Firmicutes, Deltaproteobacteria, Chloroflexi, Bacteroidetes, Nitrospirae, and Gemmatimonadetes were the core microbiome in this study. Out of all the core phyla, Alpha- and Beta-Proteobacteria were responsive to changes in soil properties and significantly predominated in SKRD compared to those of in NKRD (Supplementary Figure 1). The one-way ANOVA followed by Tukey's HSD test are shown in Supplementary Tables 4, 5.

Sixteen genera took up more than 1% of the total reads, namely Halomonas (32.10%), Aliihoeflea (5.75%), Pelagibacterium (4.89%), Pseudomonas (4.84%), Gemmata (3.56%), Nitrospira (3.13%), Lysobacter (2.65%), n class="Chemical">Nesterenkonia (2.57%), Methylocystis (2.11%), Meiothermus (1.82%), Ramlibacter (1.66%), Lactobacillus (1.56%), Methylophilus (1.34%), Bryobacter (1.30%), Bacillus (1.14%), and Burkholderia (1.04%). Core microbiome at genus level is defined as genera that existed in all the samples with at least 1% of the relative abundance within a sample. Genera Halomonas, Aliihoeflea, Pelagibacterium, Pseudomonas, Gemmata, Nitrospira, Nesterenkonia, Methylocystis, Meiothermus, Ramlibacter, Lactobacillus, Methylophilus, Bryobacter, and Bacillus were core genera in this study. Many of the core phyla had functions involved with nutrient cycling and mineral weathering in soil, which might contribute to change in soil properties.

Lefse analysis was used to identify phylum and genus that significantly responded to KRD progression (Figures 3A,B). Phylum Cyanobacteria was significantly more abundant in n class="Chemical">NKRD compared to that of in degraded soils. Genera that significantly predominated in each KRD region were determined via Lefse analysis (Figure 3A). Genera Aquicella, Rhizomicrobuim, Inquillnus, Nitrobacter, Telmatobacter, Mucilaginibacter, Anaeromyxobacter, and Rudaea were significantly abundant in NKRD areas; Roseiflexus predominated in LKRD regions; Pseudomonas, Lysobacter, Ramlibacter, Solitalea, Flavihumibacter, Microvirga, Chryseolinea, Cupriavidus, Sphingomonas, Adaeribacter, Pseudonocardia, and Roseiflexus predominated in MKRD areas; and Rhizohium, Zavarzinella, Terrimonas, Blastocatella, Altererythrobacter, Frankia, Agromyces, Actinophytocola, Senegalimassilia, and Niastella were significantly abundant in SKRD regions (Figure 3B).

Bacterial community composition in karst rocky desertification regions

The Shannon and Ace indices were calculated to assess diversity and richness of n class="Species">soil microbiome, respectively (Supplementary Figure 2). The overall change in diversity and richness were not significantly different along with KRD gradient. The diversity slightly increased; whereas richness slightly decreased in SKRD areas.

The db-RDA was used to present distribution of soil microbiome (β-diversity was calculated with unifrac-weighted matrix at phylum level) along with KRD progression (Figure 2). Result showed that n class="Species">soil microbiome in SKRD areas separated from the rest of communities along RDA1 axis, suggesting a distinguishable difference in bacterial community composition in SKRD soils compared to those of in other KRD areas. Soil microbiome in NKRD soils further separated from those of in L- and M-KRD regions along RDA2 axis, suggesting structural differences of microbiome between non-degraded and degraded soils.

Figure 2

The distance-based redundancy analysis (db-RDA) diagram shows distribution of soil bacterial communities in karst rocky desertification areas. Study sites are labeled in black No KRD (NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD). Soil properties are labeled in blue. Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Phylum is labeled in red.

The distance-based redundancy analysis (db-RDA) diagram shows distribution of soil bacterial communities in karst rocky desertification areas. Study sites are labeled in black No KRD (n class="Chemical">NKRD), latent KRD (LKRD), moderate KRD (MKRD), and severe KRD (SKRD). Soil properties are labeled in blue. Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Phylum is labeled in red.

The Linear discriminant analysis (LDA) Effect Size (LEfSe) analysis identifies phylum (A) and genus (B) that respond significantly to karst rocky desertification progression. Relative abundance is significant when P < 0.05, logarithmic LDA score ≥2.0. No KRD (n class="Chemical">NKRD).

Mantel test at 999 permutations was used to calculate impact of soil properties on bacterial community compositions across desertification gradient (Table 3). Overall, pH (r2 = 0.27, P = 0.02) was the primary factor shaping soil microbiome community in KRD regions. The MRT analysis (Figure 4) was used to interpret the relationship of relative abundance of core phylum with changes in soil properties along with KRD progression. Results indicated that soil pH was the primary factor in explaining the variation in bacterial community compositions during KRD progression. The core phyla clustered separately in soils with pH < 6.9, which only existed in n class="Chemical">NKRD regions. The MRT analysis further verified that bacterial community was mainly impacted by soil pH.

Table 3

Mantel test to determine correlation of soil properties with bacterial community composition.

Soil property	R2	F	P-value
pH	0.27	3.6	0.02*
AK	0.07	0.89	0.49
AP	0.05	0.60	0.63
TK	0.03	0.41	0.73
AN	0.03	0.32	0.80
TN	0.03	0.32	0.79
OM	0.04	0.33	0.78
TP	0.01	0.08	0.94

P < 0.05. Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Number in bold means statistically significant.

Figure 4

The multivariate regression tree (MRT) is used to reveal the relationship of relative abundance of core phyla with change in soil properties caused by KRD progression.

P < 0.05. Soil organic matter (SOM), total and available n class="Chemical">nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Number in bold means statistically significant.

Pearson Correlation analysis was used to determine correlations between soil properties and core phyla (Table 4). Soil pH and SOM had significant correlation with Acidobacteria (r2 = −0.66 and −0.69), Beta-Proteobacteria (r2 = 0.76 and 0.70), and Bacterioidetes (r2 = 0.58 and 0.54). Beta-Proteobacteria was the most sensitive phylum responding to changes in soil properties. Besides pH and SOM, other soil properties such as TN (r2 = 0.69), TK (r2 = −0.74), AP (r2 = −0.52), and AK (r2 = 0.85), also had significant correlations with Beta-Proteobacteria. Furthermore, AP (r2 = −0.54) had negative correlation with Actin class="Chemical">nobacteria.

Table 4

Pearson Correlation analysis is used to determine relationship of soil property with phylum.

phylum	pH	SOM	TN	TP	TK	AN	AP	AK
Acidobacteria	−0.26	−0.15	−0.07	0.04	−0.14	−0.10	0.29	−0.06
Alpha-Proteobacteria	−0.66**	−0.69**	−0.23	0.37	0.39	−0.39	0.48	−0.43
Planctomycetes	0.11	0.19	0.06	0.03	−0.24	0.12	−0.33	0.21
Beta-Proteobacteria	0.76**	0.70**	0.69**	−0.45	−0.74**	0.46	−0.52*	0.85***
Actinobacteria	0.25	0.26	−0.22	0.03	0.00	−0.19	−0.54*	−0.10
Firmicutes	0.40	0.36	0.39	−0.43	−0.30	0.30	−0.12	0.41
Deltaproteobacteria	−0.25	−0.32	−0.27	0.28	0.16	−0.37	0.15	−0.19
Chloroflexi	0.01	0.04	−0.21	0.07	0.17	−0.01	−0.21	−0.30
Bacteroidetes	0.58*	0.54*	0.39	−0.43	−0.31	0.39	−0.35	0.49
Nitrospirae	0.37	0.25	−0.09	0.09	0.00	−0.12	−0.68	−0.02
Gemmatimonadetes	0.37	0.41	0.19	−0.098	−0.32	0.26	−0.49	0.19

P < 0.05;

P < 0.01;

P < 0.0001. Soil organic matter (SOM), total and available nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Number in bold means statistically significant.

P < 0.0001. Soil organic matter (SOM), total and available n class="Chemical">nitrogen (TN and AN), total and available phosphorus (TP and AP), and total and available potassium (TK and AK). Number in bold means statistically significant.

Discussion

Karst rocky desertification poses as a potential threat to soil quality and ecosystem. Causes and consequences of KRD on the aboveground communities are better understood compared to those of on the underground communities, such as soil microbiome community. This study characterizes change in soil bacterial community in accordance to KRD progression, determines relationship of soil properties with bacterial community, and discusses potential activities of several bacteria in KRD regions. This study contributes to our understanding of interactions between bacterial community and KRD progression.

Change in soil properties along KRD gradient

In this study, we find that soil pH, SOM, TN, and AK significantly increase in SKRD compared to those of in n class="Chemical">NKRD areas; whereas TP and TK significantly decrease in SKRD regions. Significant increased soil pH is one of the most common indicators in karst desertification regions (Xie et al., 2015). The high soil pH in SKRD is mainly due to soil erosion and limestone dissolving. It is worth mentioning that higher SOM content alone does not indicate higher soil quality. The significantly increased SOM in SKRD areas might be resulted from excessively dissolved limestone (Zhong et al., 2006). Besides, previous studies report that the soil stability and water accessibility are low in SKRD areas, which indicate low soil quality (Peng and Wang, 2012; Qi et al., 2017). The increased TN in SKRD soils might because (1) poor vegetation coverage; and (2) predominance of nitrogen-fixing bacteria, such as Rhizohium and Frankia. Other macronutrients, TP and TK, were significant lower in SKRD than those of in NKRD areas, indicating nutrient limiting condition in SKRD soils. Interestingly, AK was significantly higher in degraded soils, which might due to the predominance and contributions of KSB, such as Pseudomonas, Bacillus, and Burkholderia (Ahmad et al., 2016). Overall, we find that soil quality in SKRD region significantly decreased compared to that of in NKRD area.

Soil pH primarily correlates with soil bacterial community

We find that soil pH is the main factor that correlates with soil bacterial community. Our results agree with numerous previous studies (Lauber et al., 2009; Rousk et al., 2010a,b; Yun et al., 2016). The overall close relationship of pH with bacterial community is due to most of bacteria having rather narrow pH optima (Charokopos et al., 2010). Slight change in pH affects more than half of bacterial activities (Fernández-Calviño and Bååth, 2010). The significant change in soil pH in this study is likely due to change in bacterial community composition. Beside direct impact, indirect effect of pH on soil bacterial community might also play an important role in shaping bacterial community. Soil pH has wide association with soil nutrient availability, catabolic activities, soil structure, and biomass activities (Wakelin et al., 2008). Lauber et al. proposed that change in soil pH altered soil fertility, structure, and vegetation community, thereby indirectly affecting soil bacterial community (Lauber et al., 2009). The broad interactions of soil pH with soil quality and habitat both directly and indirectly might explain the large impact of pH on soil bacterial community.

Many studies report responses of specific bacteria to soil pH gradient (Sait et al., 2006; Nicol et al., 2008; Jones et al., 2009). In this study, we find that soil pH significantly correlated with several core phyla (Acidobacteria, Beta-Proteobacteria, and Bacteroidetes). The correlation of soil pH with these bacteria is also reported in previous studies (Lauber et al., 2009; Yun et al., 2016). For instance, subgroups in Acidobacteria show varying preferences to soil pH. Relative abundances of subgroups 1 and 2 decrease in high pH; whereas those in subgroups 4, 6, 7, and 16 increase in high pH (Jones et al., 2009). Studies show that many bacterial activities, such as n class="Chemical">ammonia oxidization and phosphate solubilization, are pH-dependent (Nautiyal et al., 2000; Hu et al., 2013; Sharma et al., 2013). Phyla Acidobacteria, Beta-Proteobacteria, and Bacteroidetes contain many taxa carrying physiological functions that are pH-dependent, which might explain the large impact of pH has on those bacteria. Altogether, our results corroborate with previous studies, demonstrating that soil pH has a large impact on bacterial community.

Correlation of other soil properties with bacterial community

Besides soil pH, other soil properties should also be considered to fully understand the relationship of environmental factors with bacterial community. Interestingly, we find that soil pH and SOM show similar correlation with core microbiome (Acidobacteria, Beta-Proteobacteria, and Bacteroidetes), suggesting synergetic effect of pH and SOM on soil bacterial community. Indeed, a previous study reports significant correlation between soil pH and microbial decomposition rate of SOM (Yun et al., 2016). Change in soil pH affects bacterial efficiency of catabolizing n class="Chemical">cysteine, aspartic acid, lysine, and arginine (Wakelin et al., 2008). The highly alkaline soil in SKRD area results in low SOM decomposition rate, which also explain the accumulated SOM content in SKRD (Leahy and Colwell, 1990).

One member of the core microbiome, Beta-Proteobacteria, is highly responsive to change in soil properties. Six out of eight soil properties (pH, SOM, TN, TK, AP, and AK) show significant correlation with Beta-Proteobacteria. Beta-Proteobacteria contains many taxa that are highly efficient on mineral weathering in soils (Lepleux et al., 2012). Mineral weathering releases nutritive cations (e.g., n class="Chemical">phosphate, iron, and granite), which are particularly essential in nutrient limiting soils (Uroz et al., 2011). Although the result was not significant, Beta-Proteobacteria was the most abundant in M- and S-KRD soils (Figure 2). The predominance of Beta-Proteobacteria in degraded soils suggests potential nutrient releasing and cycling functions. Many taxa that belong to Bacillus (Beta-Proteobacteria) have high efficacy in mineral dissolving (Song et al., 2007; Uroz et al., 2011). The overall change in Bacillus increased along with KRD progression (Figure 4). The relatively high abundance of Bacillus further implies nutrient releasing activities in KRD area. Soil nutrient composition affects bacterial community structure (Uroz et al., 2009). For instance, various carbon sources (xylose, glucose, lactose, and mannitol) and nitrogen sources (nitrate, ammonium, and mixture of nitrate and ammonium) affect bacterial physiological activities (Nautiyal et al., 2000; Hameeda et al., 2006). The high sensitivity of Beta-Proteobacteria to change in soil properties might be explained as (1) the mineral weathering is dependent on nutrient availability in soils and (2) significant changes in soil chemical properties have large impact on Beta-Proteobacteria.

We find that phylum Cyanobacteria significantly predominate in n class="Chemical">NKRD soils. Cyanobacteria, along with lichens and mosses, compose biological soil crusts (BSCs), which preserve soil fertility in dry land (Coxson, 2002). Studies show that filamentous Cyanobacteria, such as Nostoc and Leptolyngbya, improve soil quality by providing nutrients and stabilizing soil structure (Büdel et al., 2009; Maqubela et al., 2009). Cyanobacteria provide carbon and nitrogen compounds via C and N fixing in soils (Housman et al., 2007). Cyanobacteria also produce polysaccharides that form extracellular polymeric substances (EPS) to coat and stabilize soil aggregates (Aspiras et al., 1971). The significant abundance of Cyanobacteria in NKRD regions might partially explain the better soil quality in non-degraded soils, and vice versa.

Overall, we find that soil bacterial community compositions in degraded soils are distinguishingly different from those of in non-degraded soils. Soil pH primarily correlates with soil microbiome. Other soil properties, such as SOM, Tn class="Chemical">N, TK, AP, and AK also have correlation with soil bacterial composition. Cyanobacteria, which have high potential of improving soil quality, is significantly abundant in NKRD. Many taxa with mineral weathering and nutrient cycling potentials predominate in degraded soils, suggesting that soil bacteria are important in providing nutrients in barren soils. This study fills a gap of knowledge on relationship of soil properties with bacterial community in KRD areas.

Author contributions

XW, DQ, and UD designed experiments. DQ and XZ collected and measured samples. XW analyzed data and wrote the manuscript. XW and UD edited the manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.