Heshang Cave (29°40′–30°48′ N, 108°30′–111°20′ E) is located in Changyang County, Hubei Province, China. The annual average temperature is ∼16.5°C and average annual rainfall is ∼1118 mm. Most of the precipitation (∼70%) occurs in April to September (Hu et al., 2008a; Yun et al., 2016). Vegetation consists of dense forests of subtropical tress and shrub, and corn and other vegetable crops are grown on the mountains. The cave is 250 m long and situated ∼30 m above Qingjiang River, a tributary to the Yangtze River, with a sole entrance ∼20 m in diameter (Figure 1C). The cave formed in Cambrian-aged dolomite (∼400 m thick), which is covered with a 20–40 cm thick soil (Hu et al., 2008b).

Three types of surface soils (sampling depth < 2 cm) were collected by five-point sampling method: four soil samples from the forested mountaintop (MT1/MT2/MT3/MT4), three soil samples from agricultural land (FL1/FL2/FL3) and three soil samples directly above the cave (CT1/CT2/CT3) (Figure 1C). Inside the cave, nine surface sediment samples (sampling depth < 1 cm, silt deposits) were collected from the photic zone near the cave entrance, twilight zone (i.e., the transition between the entrance and complete darkness), and the dark zone (i.e., complete darkness), with triplicate samples from each zone (Figure 1D). Nine soft weathered carbonate rock samples from the cave wall were also collected in each zone of the subterranean cave. Three drip water samples (DW1/DW2/DW3) were collected. DW3 was in the twilight zone, and DW1 and DW2 were in the dark zone (Figure 1D). Drip water samples were collected with 10 L sterile plastic bottles and sterile funnels, and the others were collected aseptically with 50 ml sterile plastic centrifuge tubes (Corning), as described previously (Yun et al., 2016). All samples were transported to the geomicrobiology laboratory in China University of Geosciences (Wuhan) under refrigeration within 24 h of collection. An aliquot of 200 ml drip water samples and a subsample of all the other samples were stored at 4°C for physicochemical property analysis. Sterile membrane filters (0.22 μm, 47 mm, Supor-200, Pall Corporation, USA) were used to filter from 3 to 10 L of the remaining water samples. The filters and remainder of other materials were stored at -80°C until DNA extraction.

pH of drip water was measured in situ with a multiparameter water quality detector (HACH, Loveland, CO, USA) (Yun et al., 2016). Solid samples were pre-frozen at -80°C for an hour and then freeze-dried (ALPHA 1-2 LD, Christ, Germany) for 48 h to remove water. An aliquot of 1 g freeze-dried solid sample was mixed with 5 ml distilled water. After shaking with an end-to-end shaker for 5 min, the mixture was centrifuged at the speed of 6,800 × g for 10 min. The supernatant pH was measured with a UB-7 pH meter (Denver Instrument) after calibrations to pH 4, 7, and 10.

Total organic carbon (TOC) content from the solid samples was analyzed with a C-S analyzer (EA 4000, Analytik Jena AG, Jena, Germany) with a 3 μg g^-1 detection limit, using high temperature ceramic technology, and standard samples AR4007 (ALPHA, USA) and AR1034 (ALPHA, USA) were used to monitor quality. TOC of drip water samples was analyzed by a TOC analyzer (Vario TOC cube, Elementar, Hanau, Germany) with a detection limit of 6 ng g^-1, using the high temperature catalytic oxidation method.

The concentrations of dissolved major anions and cations for the drip waters were measured from filtered samples. For solid samples, 2 g aliquots of freeze-dried material were mixed with 2 ml deionized water (1:1 w/v ratio) and then shaken for 10 min. After centrifugation at 2500 × g for 2 min supernatants were filtered through 0.22 μm membrane filters. The filtrates were analyzed using an ICS-600 ion chromatograph (Benedicto et al., 2014).

Total nucleic acids were extracted from 0.5 g solid sample (dry weight) using the PowerSoil DNA Kit (MoBio Laboratories, Inc., USA) following the manufacturer’s instructions. For drip water samples, the PowerWater DNA Kit (MoBio Laboratories, Inc., USA) was used to extract DNA from the membrane filters (0.22 μm, 47 mm, Supor-200, Pall Corporation, USA) according to the manufacturer’s instructions. Bacterial diversity was examined after 250-bp paired-end amplicon sequencing using the primers 520F (5′-AYTGGGYDTAAAGNG-3′) and 802R (5′-TACNVGGGTATCTAATCC-3′) (Claesson et al., 2009) on an Illumina MiSeq platform in two separate runs at Shanghai Personal Biotechnology, Co., Ltd, (Shanghai, China).

Raw sequence reads were deposited in the NCBI Sequence Read Archive under the accession number PRJNA307221¹.

Raw sequence data were quality filtered and analyzed using QIIME v 1.7.0 (Caporaso et al., 2010). Reads were processed by removing tags and primers, and the reads with an average quality score < 20 and read lengths < 150 bp were discarded. After being processed, reads were assembled by FLASH software² with the overlap between R1 and R2 reads ≥ 10 bp. Combined with mothur software v 1.31.2 (Schloss et al., 2009), chimeric sequences were identified and removed using the method of UCHIME (Edgar et al., 2011). High-quality representative sequences for each operational taxonomic units (OTUs) were assigned using UCLUST (Edgar et al., 2011) with 97% sequence identity. Taxonomic classification was carried out using Greengenes 16S rRNA database 13_8 release³ with assignment tool of Blast⁴. Sample size for each sample was rarefied prior to diversity indices calculations, including the Simpson index⁵ and Shannon index⁶, and calculations were done with mothur software by using the command of ‘summary.single.’

A heatmap of relative abundance of bacterial phyla and OTUs distribution (with relative abundance > 2%) in each sample was performed with R software (v 3.2.0)⁷. A one-way ANOVA (analysis of variance) was calculated in SPSS (version 17.0)⁸ with the significant level of 0.05 to examine significance differences among the Simpson and Shannon index values. Significant taxonomic differences between the four different habitats (e.g., overlying soils, cave sediments, weathered rocks, and drip waters) were analyzed using the least discriminant analysis (LDA) effect size (Segata et al., 2011). This method was based on the factorial Kruskal–Wallis test (α = 0.05) among classes and the pairwise Wilcoxon test (α = 0.05) between subclasses, with one-against-all strategy for multi-class analysis, to identify taxa with significant differential abundances between categories. Significant taxa were used to illustrate the difference between each sample (Bokulich et al., 2014). To demonstrate the relationship between different samples, principal coordinate analysis (PCoA) was calculated based on weighted Unifrac results in QIIME (Caporaso et al., 2010). Redundancy analysis (RDA) was performed using Canoco 5.0 to reveal the effect of environmental factors on bacterial communities, and the significant factors were chosen according to p-value (<0.05, ter Braak and Smilauer, 2002). Regression analysis and analysis of variance were used to identify the correlation between the pH and Shannon index of each sample (Anderson-Cook, 2004) in SPSS (version 17.0) software.

Table 1 summarizes the geochemical results. All solid and liquid samples were slightly alkaline, with pH ranging from 7.28 to 8.35, except four acidic mountaintop soils (Table 1). Compared with the drip waters that had low TOC values (<0.05%), TOC content was higher in overlying soils from the mountaintop and above the cave, as well as from weathered rock samples from the photic zone. NH₄⁺ was not detected in drip water or in MT and FL soils, but was detected in CT soils. NO₃^- and NO₂^- content varied in soils, with NO₃^- generally being higher than NO₂^- in most soils; except in CT soils where NO₂^- was higher than NO₃^-. Sulfate content was similar for all soil samples, ranging from 0.11 to 0.57 mM. Cave sediments and weathered rocks were chemically heterogeneous. Drip waters had higher NO₃^- and SO₄^2- concentrations compared with most of the soils, except the SO₄^2- concentration in MT4, FL2, and CT3 (Table 1).

At the phylum level, Proteobacteria and Actinobacteria were the most abundant groups in nearly all samples. For example, Proteobacteria comprised 20–40% of bacterial populations in some weathered rock samples (WA1, WA2, WA3, and WD3) and acidic surface soils (MT1, MT2, and MT3) and >60% in drip waters. Verrucomicrobia and Acidobacteria were also abundant in soils (Figure 2). Actinobacteria, Proteobacteria, and Acidobacteria were common in alkaline soils. Actinobacteria also dominated (40–60%) in all weathered rock samples and two sediment samples (SD1 and SD3). Within the Proteobacteria, Alphaproteobacteria were the most abundant in the acidic surface soils and some weathered rock samples with high TOC content (Figure 3), but their abundances decreased in the alkaline sediments and drip waters, which had the lowest TOC. Gammaproteobacteria and Betaproteobacteria were common in these low TOC samples (Figure 3).

The Simpson (Figure 4A) and Shannon indices (Figure 4B) revealed that surface soils and cave sediments had the highest bacterial diversity and drip water samples had the lowest diversity. A cluster analysis of OTU-level diversity and PCoA based on weighted Unifrac results all revealed three distinct groups (Figures 5 and 6). Group I contained bacterial communities in all the overlying soils (MT, FL, and CT) and cave sediments in photic and dark zones. Group II consisted of bacterial communities from weathered cave rocks and twilight zone sediment samples. Group III was comprised of all drip water samples, and plotted at some distance away from the other groups in PCoA space (Figure 6).

Least discriminant analysis effect size confirmed that each karst habitat had its own indicator taxa, from phylum to genus levels (Figure 7). Actinobacteria were specific to weathered rocks. Acidobacteria (18.85%), Verrucomicrobia (12.9%), and Planctomycetes (9.9%), Alphaproteobacteria (∼ 55%), and Deltaproteobacteria (∼20%) were all common in overlying soils, while Chloroflexi (9.97%), Nitrospirae (7.19%), Gemmatimonadetes (4.61%), and Firmicutes (3.53%) were common in cave sediments. Proteobacteria in drip waters were distinctive from other habitats, specifically Gammaproteobacteria and Betaproteobacteria (Figure 7). Key genera identified for each karst components were Crossiella, Euzebya, and Rubrobacter for weathered rocks, Bradyrhizobium and Acidothermus for overlying soils, Gaiella for cave sediments, and Sediminibacterium, Brevundimonas, Acidovorax, Hydrogenophaga, Polaromonas, Acinetobacter, Perlucidibaca, and Pseudomonas for drip waters (Figure 7).

At the 97% similarity level, common OTUs in weathered rocks included OTU_13, which belonged to the Pseudonocardiales (Actinobacteria), and OTU_22, which was affiliated with the Solirubrobacterales (Actinobacteria) (Figure 5). OTU_47, which belonged to the Chthoniobacterales (Verrucomicrobia), was dominant in most overlying soils (CT and MT) (Figure 5). OTU_54 of Moraxellaceae (Gammaproteobacteria) is the most abundant and unique in the drip water (Figure 5). Notably, most of the identified OTUs are putative heterotrophs according to the carbon utilization by the closest phylogenetic isolates (approximately 59.4% with relative abundances > 5‰).

To identify the environmental parameters that correlate with bacterial community variability across the Heshang Cave ecosystem, RDA was performed using 12 environmental factors (Figure 8). Among those investigated, pH (14.9% of the variability in community composition, p-value = 0.002), TOC (9.8%, p-value = 0.022) and K⁺ (7.1%, p-value = 0.03) were crucial factors that correlated with bacterial community compositions at the phylum level (Figures 8A,B). RDA analysis of proteobacterial class demonstrated 53.0% variation in total, in which pH and TOC accounted for 28.6% (Figure 8C).

Different phyla showed specific correlations with pH, TOC, and K⁺. In particular, decreased relative abundances of Acidobacteria (R²= 0.253, p-value = 0.004), Verrucomicrobia (R²= 0.5097, p-value < 0.0001) and Planctomycetes (R²= 0.3003, p-value = 0.0014) significantly correlated with an increase in pH, while higher abundances of Gemmatimonadetes and Firmicutes significantly correlated to acidic pH. Furthermore, the higher the pH, the more abundant the Gammaproteobacteria (R²= 0.167, p = 0.022) and Betaproteobacteria (R²= 0.213, p = 0.009) class representation (Figure 8C).

Relative abundances of Acidobacteria (R²= 0.228, p-value = 0.0066), Verrucomicrobia (R²= 0.1405, p-value = 0.0378), and Planctomycetes (R²= 0.2024, p-value = 0.0111) were significantly positively related to TOC content (Figure 8B). Bacterial communities from farmland soil, most of the cave sediment, and drip waters were positively correlated with pH and negatively correlated with K⁺ and TOC content in RDA space (Figure 8A). Bacterial communities in weathered rocks were plotted positively in RDA space with K⁺; Actinobacteria was the only phylum significantly correlated with K⁺ content (R²= 0.2832, p-value = 0.0021, Figures 8A,B).

Lastly, the diversity of bacterial communities in surface soils was related to pH (40.3%) (Figure 8D), and community position in RDA space indicated a strong positive correlation between the soil bacterial diversity and pH (R² = 0.7223, p-value = 0.002; Figure 9).