Microbial Diversity and Community Structure of Sulfate-Reducing and Sulfur-Oxidizing Bacteria in Sediment Cores from the East China Sea.

Sulfate-reducing bacteria (SRB) and sulfur-oxidizing bacteria (SOB) have been studied extensively in marine sediments because of their vital roles in both sulfur and carbon cycles, but the available information regarding the highly diverse SRB and SOB communities is not comprehensive. High-throughput sequencing of functional gene amplicons provides tremendous insight into the structure and functional potential of complex microbial communities. Here, we explored the community structure, diversity, and abundance of SRB and SOB simultaneously through 16S rRNA, dsrB and soxB gene high-throughput sequencing and quantitative PCR analyses of core samples from the East China Sea. Overall, high-throughput sequencing of the dsrB and soxB genes achieved almost complete coverage (>99%) and revealed the high diversity, richness, and operational taxonomic unit (OTU) numbers of the SRB and SOB communities, which suggest the existence of an active sulfur cycle in the study area. Further analysis demonstrated that rare species make vital contributions to the high richness, diversity, and OTU numbers obtained. Depth-based distributions of the dsrB, soxB, and 16S rRNA gene abundances indicated that the SRB abundance might be more sensitive to the sedimentary dynamic environment than those of total bacteria and SOB. In addition, the results of unweighted pair group method with arithmetic mean (UPGMA) clustering analysis and redundancy analysis revealed that environmental parameters, such as depth and dissolved inorganic nitrogen concentrations, and the sedimentary dynamic environment, which differed between the two sampling stations, can significantly influence the community structures of total bacteria, SRB, and SOB. This study provided further comprehensive information regarding the characteristics of SRB and SOB communities.

Introduction

Sulfur cycling, one of the key biological processes in marine sediments, is dominated by sulfate-reducing bacteria (SRB) and sulfur-oxidizing bacteria (SOB). Dissimilatory sulfate reduction, which is mediated by SRB, is considered the main process in the biomineralization of organic matter in marine sediments and might account for up to 50% of organic matter mineralization in most continental shelf sediments (Jørgensen, 1982). In addition, as much as 12–29% of the organic carbon flux on the ocean seafloor is channeled through sulfate reduction (Bowles et al., 2014). Interestingly, 80–95% of the massive amount of hydrogen sulfide formed through sulfate reduction is recycled within sediments and gradually oxidized back to sulfate (Jørgensen and Nelson, 2004). Thus, SRB and SOB control the key processes of organic matter degradation and the biogeochemical cycling of sulfur and carbon. Studying the community structure and diversity of SRB and SOB is important for revealing the roles of these bacteria in the biogeochemical cycles of carbon and sulfur and for providing insight into the biological factors driving the marine sulfur cycle.

SRB and SOB show high diversity, including both phylogenetic and metabolic diversity (Ghosh and Dam, 2009; Müller et al., 2014). To explore the environmental abundance and diversity of SRB and SOB, genes encoding key enzymes in the sulfate reduction and sulfur oxidation biochemical pathways have been used as molecular markers. For example, the dissimilatory sulfite reductase gene dsrAB, which encodes a key enzyme that catalyzes the last, essential step in the dissimilatory sulfate reduction pathway, has been frequently employed as a functional gene in studies of SRB in various environments (Foti et al., 2007; Pester et al., 2012; He et al., 2015). The known SOB have been demonstrated to use different enzymes, pathways, and electron transport and energy conservation mechanisms for the oxidation of sulfide. The sulfur-oxidizing (Sox) (Meyer et al., 2007), reverse dissimilatory sulfite reductase (Loy et al., 2009), adenosine-5-phosphosulfate reductase (Meyer and Kuever, 2007), and sulfide quinone oxidoreductase enzyme systems (Pham et al., 2008) play vital roles in sulfide oxidation, and among these, the Sox multi-enzyme system is considered a fundamental and primordial molecular mechanism for sulfur oxidation (Ghosh and Dam, 2009) that is widespread among the known SOB (Meyer et al., 2007). Moreover, the soxB gene, which encodes the SoxB subunit of the Sox enzyme system, has been widely employed to characterize the abundance and diversity of SOB in various environments (Kojima et al., 2014; Thomas et al., 2014; Tourna et al., 2014).

The diversity and community structure of SRB and SOB have, to date, been studied using denaturing gradient gel electrophoresis (DGGE) (Varon-Lopez et al., 2013; Yang et al., 2015), restriction fragment length polymorphism (Luo et al., 2011; Reed and Martiny, 2013), and clone library approaches (Yang et al., 2013; Purcell et al., 2014). However, obtaining comprehensive information about the highly diverse SRB and SOB communities, particularly SRB and SOB from the prime habitats of marine sediments, using these methods is challenging (Aoki et al., 2015; Zhang et al., 2016). High-throughput sequencing, which constitutes a powerful approach for achieving complete coverage of microbial communities, has been recently applied in the analysis of microbial community diversity and composition. For instance, high-throughput sequencing of the 16S rRNA gene revealed the presence of a number of novel bacterial groups (Zhu et al., 2013; Mahmoudi et al., 2015) and a high diversity of bacteria (Wang et al., 2015). Moreover, the community structure and diversity of functional microbes, particularly nitrogen-cycling microbes, such as ammonia-oxidizing bacteria, ammonia-oxidizing archaea and denitrifying microbes, have been investigated through high-throughput sequencing analyses of functional genes (Tago et al., 2014; Saarenheimo et al., 2015). However, only a few studies have evaluated the community structure and diversity of SRB and SOB using high-throughput sequencing (Meyer et al., 2016; Cui et al., 2017).

The East China Sea (ECS) is the largest marginal sea of the Northwest Pacific Ocean, with a vast continental area of 0.5 × 1012 m2. The Changjiang River (Yangtze River), which transports great quantities of freshwater, nutrients, organic matter, and other chemical elements into the ECS, has a marked effect on the sea (Milliman and Meade, 1983; Beardsley et al., 1985). The ECS is also influenced by the warm and oligotrophic Taiwan Warm Currents in the south and the Kuroshio Current (from the Western Equatorial Pacific) in the east (Jiao et al., 2005). More interestingly, the ECS has also experienced serious environmental problems, including eutrophication, harmful algal blooms, and chemical pollution (Shen et al., 2011; Cheng et al., 2012; Yu et al., 2013). Because of the strong influence of its riverine inputs, currents and environmental problems, the ECS exhibits complex physical and chemical conditions and has become an ideal study area for ecological investigations of the temporal and spatial dynamics of biota (Dang et al., 2008; Xiong et al., 2014; Zhang et al., 2014). Increasing research efforts have focused on microbial ecology in the ECS, particularly the microbial community structure and diversity of organisms associated with the sulfur cycle, such as SRB and SOB (Nie et al., 2009; Wu et al., 2009; He et al., 2015); however, because of the limitations of the applied methods (DGGE or clone libraries), only a few dominant species have been identified. In the present study, the 16S rRNA gene and two functional genes (dsrB and soxB) were combined in a high-throughput sequencing analysis, and this combination allowed for an in-depth analysis of the community structure and diversity of SRB and SOB in sediment cores from the ECS. In addition, quantitative PCR (qPCR) was performed to reveal the vertical distribution of SRB and SOB populations and provide insight into the SRB and SOB community characteristics.

Materials and Methods

Site and Sampling Description

Two core samples from ECS sediments were collected from station S31 (length: 21 cm) and station S33 (length: 50 cm) during a research cruise in July 2011 on the R/V Run-Jiang (Supplementary Figure S1). A QuAAtro nutrient auto analyzer (Seal Analytical Ltd., United Kingdom) was used to measure the pore-water dissolved N, P, and Si concentrations of the following: nitrate (NO3-N), nitrite (NO2-N), ammonium (NH4-N), phosphate (PO4-P), and silicate (SiO3-Si) (Supplementary Table S1). The sampling methods as well as the methods used for the determination of environmental parameters [pH and the concentrations of Fe(II), Mn(II), sulfate and excess 210Pb] were described in detail by He et al. (2015) (Supplementary Table S1 and Figure S2). After collection, the samples from stations S31 and S33 were divided into three depth sections (0–4, 8–12, and 16–20 cm) and five depth sections (0–4, 8–12, 16–20, 32–36, and 46–50 cm), respectively, and subsequently frozen at -80°C for nucleic acid extraction.

DNA Extraction

Genomic DNA was extracted from each sediment sample un class="Chemical">sing the PowerSoil DNA Isolation Kit (Mo Bio Laboratories Inc., Carlsbad, CA, United States) according to the manufacturer’s recommended protocol. The extracted DNA was stored at -80°C prior to further analyses.

High-Throughput Sequencing and Data Processing

Illumina HiSeq2500 (PE250) Sequencing and Analysis

The community structure, richness, and diversity of total bacteria and SRB were studied through Illumina HiSeq2500 (PE250) sequencing based on the 16S rRNA and dsrB genes. Fragments of the 16S rRNA gene (∼466 bp, V3–V4 region) and dsrB gene (∼350 bp) were amplified using the primer pairs 341F/806R and DSRp2060F/DSR4R, respectively (Supplementary Table S2); these primer pairs have been widely employed in previous studies of bacteria from various environments (Geets et al., 2006; Michelsen et al., 2014). PCR amplification was performed in a 30-μL reaction volume containing 15 μL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM forward and reverse primers, and approximately 10 ng of template DNA. The amplification reactions were performed in a thermal cycler (Bio-Rad T100, United States) and consisted of an initial denaturation step at 98°C for 1 min followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, elongation at 72°C for 30 s, and a final step at 72°C for 5 min. The PCR products were analyzed via 2% agarose gel electrophoresis to assess the quality and size of the resultant amplicons. The PCR products were then mixed at equidensity ratios and purified using the GeneJET Gel Extraction Kit (Thermo Fisher Scientific). 16S rRNA and dsrB gene libraries were subsequently prepared using the NEB Next®UltraTM DNA Library Prep Kit (New England Biolabs), and the libraries were sequenced on the Illumina HiSeq2500 (PE250) platform following the manufacturer’s recommendations.

Paired-end reads (16S rRNA and dsrB) were assigned to the samples based on their unique barcode, truncated by removing the barcode and primer sequence, and then merged using FLASH (Version 1.2.7) (Magoè and Salzberg, 2011). The merged reads were subsequently subjected to quality-based filtering with QIIME (Caporaso et al., 2010). Briefly, the reads were truncated at any site containing more than three sequential bases with a Phred quality score (Q) below 20 and any read containing ambiguous base calls and reads with less than 75% (with respect to the total read length) consecutive high-quality base calls (Wang et al., 2016) were discarded.

Pyrosequencing and Sequence Analysis

The soxB gene (∼750 bp) was amplified with the soxB693F/soxB1446B primer set, which was previously shown to yield the most successful and reliable amplification results (Supplementary Table S2) (Meyer et al., 2007). The 25-μL amplification reaction for the soxB gene included 4 μL of 5 × Q5 reaction buffer (New England Biolabs), 2 μL of 2.5 mM deoxynucleoside triphosphate (dNTP) mix, 5 μL of 5 × Q5 High GC Enhancer, 1 μL of each primer (10 μM), 0.25 μL of Q5 High-Fidelity DNA Polymerase (New England Biolabs), 2.5 μL of template DNA, and 8.25 μL of double-distilled H2O. The reactions were maintained at 98°C for 5 min for DNA denaturation; they were then subjected to 32 cycles of 98°C for 30 s, 55°C for 40 s, and 72°C for 1 min and then to a final extension at 72°C for 7 min to ensure complete amplification. The obtained PCR product was cut from a 1.5% agarose gel and purified using the AxyPrep DNA Gel Extraction Kit (Axygen, AP-GX-250). Pyrosequencing of the soxB gene was subsequently conducted using the 454 FLX+ platform according to the manufacturer’s recommended protocol.

The raw data were processed un class="Chemical">sing the QIIME pipeline (Caporaso et al., 2010). For quality filtering, low-quality reads shorter than 150 bp, reads with an average quality score below 19, sequences with ambiguous base calls and sequences with homopolymer runs exceeding 6 bp were removed (El-Chakhtoura et al., 2015). Moreover, the barcode and primer regions were trimmed from the sequences.

Operational Taxonomic Unit Clusters and Taxonomic Assignment

After quality control, UPARSE (Edgar, 2013) was employed to cluster all of the clean reads into operational taxonomic units (OTUs) using 97 and 90% similarity cutoffs (Petri et al., 2001; Müller et al., 2014) for the 16S rRNA gene and functional genes (dsrB and soxB), respectively. In addition, the most abundant sequence from each OTU was selected as the representative sequence. Taxonomic information for each representative 16S rRNA and functional gene sequence was obtained by searching the Greengenes and the NCBI (National Center for Biotechnology Information) databases, respectively. Reads that did not match any sequences in the database were clustered into the unclassified group.

Quantification of Gene Copies in Sediments

The total bacterial 16S rRNA gene, the SRB dsrB gene, and the SOB soxB gene in the sediment samples were amplified using the primers 341F/518R, DSRp2060F/DSR4R, and soxB693F/soxB1164BK145, respectively (Supplementary Table S2); these primers were demonstrated to work effectively in previous studies (Geets et al., 2006; Dang et al., 2010; Krishnani et al., 2010). The method used for quantification of the dsrB and 16S rRNA genes in the investigated samples has been described previously (He et al., 2015). All qPCR assays targeting the soxB gene were performed in triplicate using an ABI PRISM® 7500 Sequence Detection System (Applied Biosystems, Foster City, CA, United States). Each 20-μL qPCR mixture contained the following reagents: 10 μL of FastStart Universal SYBR Green Master Mix (Rox) (Roche Diagnostics, Mannheim, Germany), 0.2 μg μL-1 of bovine serum albumin, 0.4 μM of each primer, and 2.0 μL of template DNA. The qPCR amplification conditions were as follows: 95°C for 10 min followed by 40 cycles of 15 s at 95°C, 45 s at 55°C, and 1 min at 72°C. To evaluate the specificity of the qPCR amplification, the PCR products were sequenced by the Beijing Genomics Institute. In addition, the specificity of the amplification products was verified by melting curve analysis and visualized in agarose gels. Melting curves were obtained at 60–95°C; a read was obtained after every 1°C increase in temperature, and the temperature was maintained for 1 s between reads. The resultant qPCR data were analyzed using ABI PRISM 7500 SDS software.

Plasmids containing the target gene fragments were extracted from Escherichia coli hosts using a FastPlasmid Mini Kit (CWBIO, Beijing, China) and quantified using a Picodrop microliter spectrophotometer (Picodrop, Saffron Walden, Essex, United Kingdom). Standard curves for the qPCR assays were obtained with serial 10-fold dilutions of reference plasmids. Standard curves with efficiencies ranging from 90 to 110% and a corresponding R-value greater than 0.99 were considered credible.

Statistical Analyses

After randomly reducing the number of reads to the lowest number of reads in any individual sample, the richness indices (Chao 1 estimates), diversity indices (Shannon index), and Good’s coverage were obtained with QIIME (Caporaso et al., 2010). Furthermore, unweighted pair group method with arithmetic mean (UPGMA) clustering based on a weighted UniFrac distance matrix was conducted using QIIME to examine the differences in the bacterial communities between the sediment samples. The relationships between community composition and sediment characteristics were examined using BIOENV, detrended correspondence analysis (DCA), and redundancy analysis (RDA) with the R package vegan (Oksanen et al., 2011). The DCA analysis indicated that a linear-model-based RDA analysis was more suitable than unimodal CCA for our data (Hernández-Landa et al., 2015). BIOENV provides the subset of environmental variables that best explain the community variation among the sampling stations using Spearman rank correlation, and RDA determines the percentage of the community composition variation explained by this subset of environmental variables. The significance of the environmental variables was tested by a Monte Carlo permutation test (999 unrestricted permutations, p < 0.05). A Pearson’s correlation analysis between microbial abundance and environmental parameters was performed using SPSS statistics software.

Nucleotide Sequence Accession Numbers

The sequence data generated in this study were depon class="Chemical">sited in the NCBI Short Read Archive database under accesn class="Chemical">sion numbers SRP077048 (16S rRNA), SRP077077 (dsrB), and SRP077086 (soxB).

Results

Bacterial 16S rRNA Gene Analysis

A total of 371,534 high-quality 16S rRNA gene sequences with an average length of 418 bp were obtained from an Illumina HiSeq2500 (PE250) sequencing analysis of the eight sediment samples. The high-quality sequences were clustered into OTUs with 97% sequence similarity. Accordingly, bacterial richness (Chao 1 estimator) and diversity (Shannon index) estimates were calculated for each sample and are shown in Table . Good’s coverage estimate for each sample ranged from 98.2 to 99.2%, indicating that the sampling was sufficient to cover almost all bacterial communities. A total of 1,607–2,196 OTUs with 97% similarity were identified in the eight sediment samples. The samples from S33 exhibited more OTUs than those from S31, and the surface sediment (0–4 cm) samples presented the lowest number of OTUs at each station. The diversity indices (Shannon indices) and species richness indices (Chao 1) indicated that the samples from S33 displayed greater richness and diversity than those from S31. These results revealed higher microbial diversity at S33. Furthermore, the diversity indices did not change markedly with increasing depth. In addition, the lowest and the highest richness indices were obtained for the 0–4 and 32–36 cm sections, respectively.

Vertical distribution of sediment bacterial community divern class="Chemical">sity and richness estimators based on the n class="Chemical">16S rRNA gene.

In total, 52 different bacterial phyla were detected across all sediment samples, and the dominant phyla (exhibiting a relative abundance >1% in at least one sample), which accounted for more than 94.57% of the total sequences in each sample, are shown in Figure . Proteobacteria was the most abundant phylum in all samples and accounted for 36.9–62.4% of all sequences. After the relative abundance of Proteobacteria reached its peak value at a depth of 8–12 cm, it exhibited a tendency to significantly decline with increasing depth at each station. Within Proteobacteria, the majority of sequences were assigned to Deltaproteobacteria (25.8–70.0%), Gammaproteobacteria (12.7–46.5%), and Alphaproteobacteria (11.1–34.3%) (Figure ). Furthermore, Deltaproteobacteria and Alphaproteobacteria were more abundant at S33, whereas Gammaproteobacteria was more abundant at S31.

n class="Chemical">16S rRNA gene-based relative abundance of dominant bacterial groups (relative abundance >1%) observed in core samples collected from the East China n class="Gene">Sea.

n class="Chemical">16S rRNA gene-based relative abundance of different classes of Proteobacteria in core samples collected from the East China n class="Gene">Sea.

In the shallow sediment layers (0–20 cm) from each station, Actinobacteria and Acidobacteria were the second and third most dominant communities, accounting for 5.8–27.6 and 6.5–9.4% of all bacteria amplicons, respectively. Furthermore, the relative abundance of Actinobacteria decreased with increasing depth at each station, whereas the relative abundance of Acidobacteria increased with increasing depth at S31 and decreased with increasing depth at S33. In the deeper sediment layers (32–36 and 46–50 cm) from S33, Chloroflexi and Crenarchaeota were the second and third most dominant communities, accounting for 13.1 and 11.4% of bacterial amplicons, respectively.

Above all, the bacterial communities at the two investigated stations showed similarities in their dominant groups of phyla, particularly in shallow sediment layers. However, certain variations, such as differences in class groups, were observed between the two stations. The similarities and dissimilarities between the bacterial communities in the different sediment samples were further quantified through UPGMA clustering analysis based on the weighted UniFrac distance metric (Figure ). Overall, the bacterial community structures in the shallow sediment layers (0–20 cm) from S31 and S33 were similar, and those in the deeper sediment layers (32–36 and 46–50 cm) from S33 clustered separately. These findings, which showed high jackknife support, suggest that sampling depth plays a relatively important role in bacterial community structure. However, analysis of the cluster containing the shallow sediment samples revealed that the samples from S31 formed a separate group and clustered away from the samples from S33, revealing the existence of differences in bacterial community structure between the two sampling sites.

Tree representing the results of the UPGMA hierarchical clustering of the weighted UniFrac distance matrix for the 16S rRNA (A), dsrB (B), and soxB (C) genes in sediments from the ECS. The scale bar indicates the distance between clusters in UniFrac units. The red branch represents 75–100% level of jackknife support.

Functional Gene Analysis

Functional Gene Diversity

The total numbers of high-quality sequences from the dsrB and soxB genes obtained from all samples were 615,062 and 107,688, respectively. Using a similarity cutoff of 90%, these sequences were classified into 829 OTUs for the dsrB gene and 380 OTUs for the soxB gene. The corresponding species richness (Chao 1 estimator) and diversity estimates (Shannon index) were calculated for each sample and are listed in Tables . Good’s coverage values of the functional genes were all above 99%, indicating that the obtained sequences could adequately reflect the diversity of SRB and SOB in the sediments. Overall, the soxB gene diversity was substantially higher than that of the dsrB gene, whereas the dsrB gene richness was markedly higher than that of the soxB gene.

n class="Chemical">Similarity-based OTUs and species richness and divern class="Chemical">sity estimates based on the dsrB gene.

n class="Chemical">Similarity-based OTUs and species richness and divern class="Chemical">sity estimates based on the soxB gene.

The diversity indices of the dn class="Chemical">srB gene did not change markedly with increasing depth at S31 and increased with increasing depth at S33. Moreover, the diversity indices of the dsrB gene in the shallow sediment layers (0–20 cm) from S31 were greater than those from S33 and lower than those obtained for the deeper sediment layers (32–36 and 40–46 cm) from S33. The richness of the dsrB gene was lowest in the 8–12 cm sediment depth range at station S31 and reached its peak value in the depth range of 32–36 cm at station S33. For the soxB gene, the diversity and richness indices decreased with increasing depth at S31, whereas at S33, the diversity and richness indices decreased with increasing depth until suddenly increasing and peaking in the depth range of 46–50 cm.

Functional Gene Community Composition

dsrB gene

The community structure of the dsrB gene in the samples was analyzed based on the 109 dominant (showing >1% abundance in at least one sample) and core (common to all samples) OTUs (representing 85.3–98.2% of the total sequences). Overall, 103 OTUs, accounting for 81.1–94.7% of the total dsrB sequences, were affiliated with Deltaproteobacteria, and six OTUs, representing 0.45–6.6% of the total dsrB sequences, were classified as belonging to Firmicutes. As shown in Figure , among the dominant Deltaproteobacteria, the majority of sequences (15 OTUs with relative abundances ranging from 40.91 to 73.36%) could not be assigned at the class or family levels. The 42 OTUs belonging to the Desulfobacteraceae family (ranging from 9.60 to 34.05%) represented a significant fraction and tended to increase with increasing depth at each station. A total of 42 OTUs with relative abundances ranging from 4.23 to 7.49% were affiliated with Syntrophaceae. The Desulfobulbaceae family (four OTUs) reached an abundance of 4.06% on average at S31 but only accounted for 0.40% of the sequences at S33.

Community structure of the dn class="Chemical">srB gene in sediment samples at the family level.

The relative abundance of all dominant OTUs (total of 28 OTUs) from each sample exceeded 80%, indicating that the community composition of the dsrB gene in sediments was well reflected by the dominant OTUs. Thus, the distribution of dominant OTUs was analyzed to better understand the composition and structure of the dsrB gene in the samples (Figure ). OTU1, showing relative abundances ranging from 38.8 to 63.8%, was the predominant SRB group in all sediment layers. Moreover, the relative abundance of OTU1 decreased with increasing depth at S33, whereas at S31, no significant difference was observed between the depth ranges of 0–4, 16–20, and 8–12 cm. In addition, the relative abundance of OTU1 in the shallow sediment layers (0–20 cm) from S33 was greater than that in the samples from S31. OTU2 reached an average abundance of 3.98% at S31 but only accounted for 0.07–0.5% of the sequences at S33. A similar trend of differences was observed for OTU5, OTU6, OTU7, OTU8, and OTU9, whereas OTU3, OTU4, OTU10, OTU11, and OTU12 were more abundant at S33 than at S31.

Relative abundance of dominant OTUs (relative abundance >1%) of the dn class="Chemical">srB gene observed in core samples collected from the East China n class="Gene">Sea.

Taken together, the results regarding the distribution of dominant OTUs at S31 and S33 indicated certain differences. Interestingly, similar results were obtained via UPGMA clustering analysis (Figure ). The UPGMA analysis showed that the SRB communities at S31 and S33 clustered separately from each other, revealing the dissimilarity in the SRB community structure between the two sampling stations and among different sediment depths.

soxB gene

Taxonomic classification revealed that Proteobacteria accounted for the most abundant clades at the phylum level, representing 46.91–77.17% (131 OTUs) of the sequences. Chlorobi and Spirochaetes were the second and third most dominant phyla, accounting for 5.32 (two OTUs) and 2.31% (eight OTUs) of the sequences, respectively (Figure ). In addition, 238 OTUs representing 20.05–45.82% of the sequences were n class="Species">unclassified (Figure ). The dominant Proteobacteria phylum was classified into three classes: Alphaproteobacteria (varying from 31.37 to 53.00%), Betaproteobacteria (varying from 1.00 to 14.64%), and Gammaproteobacteria (varying from 4.28 to 12.50%).

Relative abundance of the n class="Gene">soxB gene observed in core samples collected from the East China n class="Gene">Sea.

The composition and structure of the SOB community in the various samples were compared at the family level (Figure ). Bradyrhizobiaceae was the dominant family in all samples and was more abundant (34.72%) in S31 than in S33 (ranging from 15.79 to 29.86%). Four other minor families (Rhodobacteraceae, Hyphomicrobiaceae, Rhodospirillaceae, and an unclassified group) within Alphaproteobacteria were more abundant in S33 than in S31. Among the Betaproteobacteria, Burkholderiaceae, with a relative abundance of 4.31% (on average), was the dominant family in S31, whereas an unclassified group with a relative abundance of 5.91% (on average) was the dominant family in S33. Furthermore, within Gammaproteobacteria, Ectothiorhodospiraceae was the dominant family in almost all samples from S33, but its abundance was only 0.11–0.71% at S31. The UPGMA clustering analysis (Figure ) based on the weighted UniFrac distance metric showed that the SOB communities from S31 and S33 clustered separately from each other, further demonstrating the differences in the SOB community structure between the two sampling stations and among different sediment depths.

Relative abundance of the n class="Gene">soxB gene in the core samples at the family level.

Influence of Environmental Variables on Total Bacterial, SRB, and SOB Communities

BIOENV identified that the environmental factors most strongly correlated with the total bacteria, SRB, and SOB communities were depth, pH, dissolved inorganic nitrogen (DIN), phosphate (PO4-P), and silicate (SiO3-Si) (correlation >0.6). Together, these variables explained 76.58, 90.10, and 88.25% of the variation in the bacteria, SRB, and SOB communities, respectively, as determined through RDA (Figure ). The environmental variables that were found to contribute significantly to the microbial community–environment relationship (900 Monte Carlo permutations) were depth and DIN for all of the examined genes (16S rRNA, dsrB and soxB gene).

Redundancy analyn class="Chemical">sis (RDA) biplot of environmental parameters and the n class="Chemical">16S rRNA and functional genes (dsrB and soxB). (A) 16S rRNA; (B) dsrB; (C) soxB.

Abundance of Total Bacterial 16S rRNA and Functional Genes

The abundance of total bacterial 16S rRNA and functional genes (dsrB and soxB) in the 21- and 50-cm-long core samples from S31 and S33, respectively, were determined via qPCR (Figure ). In the core sample from S31, the vertical abundance profile of the dsrB, soxB, and 16S rRNA genes showed marked fluctuations, and the range of fluctuation in the abundances of the 16S rRNA and soxB genes was smaller than that in the abundance of the dsrB gene. However, analysis of the core sample from S33 showed that the abundances of the dsrB, soxB, and 16S rRNA genes gradually decreased with increasing depth.

Vertical profiles of the dn class="Chemical">srB, n class="Gene">soxB, and total bacterial 16S rRNA gene copy numbers in the core samples from S31 and S33.

Discussion

In this study, a comparative analysis of the sulfur cycle-related microbial communities (SRB and SOB) in sediment cores was performed using an approach combining 16S rRNA, dsrB, and soxB gene high-throughput sequencing with qPCR analyses. The combination of these molecular techniques provided a detailed description of the SRB and SOB communities. The traditional approaches (e.g., DGGE and clone library), which have been widely applied for the detection of SRB and SOB in various environments, might notably underestimate their diversity. High-throughput sequencing based on functional genes (dsrB and soxB) can result in a greater number of sequences and almost complete coverage (Good’s coverage values were all greater than 99%), providing more comprehensive information regarding the diversity and community composition of the SRB and SOB communities. Thus, the diversity, richness and OTU numbers observed in our high-throughput sequencing analyses were notably greater than those obtained in previous studies that employed other molecular methods (Jiang et al., 2009; Yang et al., 2013; Yousuf et al., 2014; He et al., 2015). Further analysis indicated that rare species make important contributions to the greater richness and diversity of the SRB and SOB communities. For instance, as shown in Table , the number of dominant OTUs in each sample was similar, ranging from 8 to 13 and representing 3.08% of OTUs (on average); however, these OTUs comprised more than 77% of the sequence abundance for the dsrB gene. Although the rare OTUs (relative abundance <0.01%) (Galand et al., 2009) accounted for 54.03% of OTUs (on average), they only comprised 0.84% (on average) of the sequence abundance. Similar results were obtained for the soxB gene (Table ).

Rare and dominant OTUs based on the dn class="Chemical">srB gene.

Rare and dominant OTUs based on the n class="Gene">soxB gene.

The 16S rRNA sequencing analysis identified 235 and 191 OTUs as belonging to potential SRB and SOB, respectively, based on family species classification. Furthermore, these SRB OTUs were affiliated with 10 families, accounting for 7.63–19.62% (14.12% on average) of the total 16S rRNA gene sequences (Table ), and the SOB OTUs were affiliated with 12 families, accounting for 5.99–18.78% (12.53% on average) of the total 16S rRNA sequences (Table ). This finding is consistent with the results of a previous study in this region that found that numerous sequences were affiliated with potential SRB and SOB (Ye et al., 2016). Overall, these results confirm the high abundance of potential SOB and SRB groups in the ECS sediments and suggest the existence of an active sulfur cycle in this area, and these sulfur cycle-related bacterial communities have considerable potential for further exploration.

n class="Chemical">SRB community compon class="Chemical">sition based on taxonomic information of the dsrB and 16S rRNA genes.

n class="Chemical">SOB community compon class="Chemical">sition based on taxonomic information of the soxB and 16S rRNA genes.

The dsrB gene sequencing analysis performed in this study identified two phyla (Proteobacteria and Firmicutes), and Deltaproteobacteria (within Proteobacteria) represented the dominant SRB. Our data were consistent with previous observations obtained from various environments, such as marsh sediments (Jemaneh et al., 2013), the Great Salt Lake (Kjeldsen et al., 2007), mangrove sediments (Varon-Lopez et al., 2013), wastewater treatment plants (Biswas et al., 2014), and marine and estuary sediments (Jiang et al., 2009; Wu et al., 2009; He et al., 2015), suggesting that Deltaproteobacteria have a wide adaptation range and play major roles in global sulfur cycling. The soxB gene sequencing results revealed that the SOB community was dominated by Alphaproteobacteria and unclassified groups, followed by Gammaproteobacteria and Betaproteobacteria. For comparison, Gammaproteobacteria and Epsilonproteobacteria are considered dominant Sox organisms in marine sediments (Sievert et al., 2008; Lenk et al., 2011; Akerman et al., 2013; Dyksma et al., 2016); however, no Epsilonproteobacteria groups were observed in the present study. This discrepancy could be due to primer bias; the primers used for soxB gene amplification amplify the soxB genes of Epsilonproteobacteria and Chloroflexi very poorly (Meyer and Kuever, 2007). In addition, Bradyrhizobiaceae within Alphaproteobacteria and Burkholderiaceae within Betaproteobacteria were the dominant families in our study area; however, previous research revealed that these families might adapt to oligotrophic sulfur environments, such as soil (Jung et al., 2005; Masuda et al., 2010) and rhizosphere soils (Anandham et al., 2008; Masuda et al., 2016). Therefore, further study is required to determine whether the families described here are the result of contamination or actually exist. More interestingly, most of the domain and core OTUs obtained through dsrB and soxB gene high-throughput sequencing were difficult to clearly assign at the family or genus level, partially due to the limited number of reference sequences. For instance, 15 OTUs for the dsrB gene, with a relative abundance ranging from 40.91 to 73.36%, and 238 OTUs for the soxB gene, with a relative abundance ranging from 20.05 to 45.82%, could not be assigned at the family or genus level. Overall, these findings suggest that the ECS environment contains a high abundance of yet unknown SRB and SOB lineages, which should be studied in more detail. Indeed, this phenomenon is consistent with the following recent findings: 62% of dsrB sequences in surface sediments from the ECS exhibit low similarity with previously cultured SRB (Zhang et al., 2016); 57.6% of dsrAB clones from polluted harbor sediments belong to unclassified groups of SRB (Zhang et al., 2008); and 60% of sequences obtained through a dsrA gene analysis show no clear affiliation with a known SRB (Quillet et al., 2012). Furthermore, similar results have been reported by Yousuf et al. (2014) for soxB. In addition, we noticed that these unclassified OTUs exhibit high similarity with environmental dsrB/soxB sequences, as demonstrated through BLAST identity searches. For example, most of the unclassified dsrB OTUs displayed high similarities (>90%) with uncultured SRB recovered from Japanese fish farm sediments (Kondo et al., 2012) and ECS sediments (He et al., 2015; Zhang et al., 2016). Similarly, most of the unclassified soxB sequences exhibited high similarity (>80%) with sequences recovered from coastal soil ecosystems (Yousuf et al., 2014). These observations suggest that these unclassified sequences are widely distributed in marine and terrestrial environments and play a vital role in the biogeochemical cycling of sulfur and carbon.

A comparison of the SRB community composition identified based on the 16S rRNA and dsrB genes showed that 16S rRNA gene sequencing detected 10 SRB families, whereas dsrB gene sequencing only detected four SRB families (Tables ). Moreover, it was difficult to obtain more detailed taxonomic information for the SRB community based on the 16S rRNA and dsrB genes. Overall, the dsrB gene proved to be less efficient for the proper annotation of taxonomic information for SRB groups. This phenomenon might be attributed to the limited number of reference sequences used to annotate taxonomic information for the dsrB gene. In contrast, the results obtained from the high-throughput sequencing of the soxB gene revealed 15 SOB families, whereas the high-throughput sequencing of the 16S rRNA gene revealed 12 SOB families. In addition, the soxB gene can provide more detailed taxonomic information (at the genus or even the species level) for the SOB community (data not shown) than the 16S rRNA gene. Taken together, the results indicate that the soxB gene might better define SOB groups.

The 16S rRNA, dsrB and soxB gene high-throughput sequencing data obtained in the present study revealed a high diversity of SRB and SOB in sediments from the ECS, and the classification of SRB and SOB communities based on 16S rRNA gene and functional gene was analyzed and compared according to the results of a previous study (Meyer et al., 2016). Overall, the present study provided more comprehensive information regarding the SRB and SOB community characteristics and suggested the existence of an active sulfur cycle in the study area. Significantly, the sulfate amount never decreases below approximately 24 mM in the core samples, which appears to contradict the conclusion that SRB and SOB communities are responsible for active sulfur cycling. In fact, this phenomenon will be observed based on the fact that the sulfate consumed by sulfate reduction is recycled within sediments by oxidation of reduced inorganic sulfur compounds to sulfate. As previously reported, 80–95% of the massive amount of hydrogen sulfide formed through sulfate reduction is recycled within sediments and gradually oxidized back to sulfate (Jørgensen et al., 1990; Jørgensen and Nelson, 2004). Jørgensen et al. (1990) reported that sulfate reduction rates increased in the top 10 cm of core in the Belt Sea; however, the concentration of sulfate varied little (approximately 25 mM). Similarly, Leloup et al. (2009) revealed that substantial sulfate reduction rates were measured in the top 25 cm of core, and the sulfate concentration profile was approximately 20 mM. However, the diversity of sulfur cycle-related bacteria and the potential sulfur cycle activity might still be overestimated or underestimated because certain named SRB/SOB might did not actually take part in the sulfur cycle in sediments, such as Pelotomaculum and Sporotomaculum, which are not able to grow with sulfite and/or sulfate as electron acceptors even though they are closely affiliated with SRB (Brauman et al., 1998; Imachi et al., 2006). Instead, certain unnamed SRB/SOB, such as archaeal anaerobic methanotrophs, might be involved in the sulfur cycle in sediments (Treude et al., 2014). Fortunately, these groups are only found in special habitats and account for a small portion of the identified communities (Brauman et al., 1998; Imachi et al., 2006; Miyashita et al., 2009; Treude et al., 2014). Thus, the use of 16S rRNA and functional genes (dsrB and soxB) as proxies for determining the microbial groups that play active roles in sulfur cycling and for studying the diversity of SOB and SRB is a powerful approach.

The sediment environments of the ECS are highly heterogeneous (Dang et al., 2008; Fang et al., 2013); this characteristic has a strong influence on the spatial heterogeneity of the sediment microbial communities in the ECS, as observed in previous studies (Dang et al., 2008; Feng et al., 2009; Chen et al., 2014; Ye et al., 2016). Based on the sediment grain size combined with the overlying water masses, the ECS can be divided into three domains: the inner shelf mud area, the outer shelf sand area, and the slope plus Okinawa Trough mud area (Liu et al., 2006; Xu et al., 2016). S31 is located within the inner shelf mud deposits along the Zhejiang coast, whereas S33 is located at the periphery of the Zhejiang coastal mud area (Supplementary Figure S1). Thus, most of the environmental parameters showed substantial variation among the different samples and among the various sediment depths (Supplementary Table S1 and Figure S2). The concentration of DIN, for instance, ranged from 91.44 to 449.21 μM in station S31 and from 18.81 to 120.73 μM in station S33 (Supplementary Table S2). 210Pb has been proven to serve as an effective proxy for sedimentary dynamics processes, such as erosion–transportation–deposition, boundary scavenging, and resuspension process, in coastal and shelf regions (Huh and Su, 1999; Su and Huh, 2002). The excess 210Pb activity varied markedly between the two sediments (Supplementary Figure S2). At station S31, 210Pbex activity was stable from the surface sediment to a depth of 8 cm; however, at a depth below 12 cm, 210Pbex activity increased gradually and displayed characteristics of reverse accumulation, indicating that this sedimentary layer has experienced sedimentary dynamic events. At station S33, 210Pbex activity remained stable until a depth of 5 cm and then gradually declined. This finding indicates that the sedimentary dynamic environment at S33 was more stable than that at S31, as observed in our previous study (He et al., 2015). Our results showed that the abundances of 16S rRNA, dsrB and soxB genes show substantial fluctuations at S31 (Figure ), and none of the gene abundance values were significantly correlated with sediment depth. In contrast, all of the genes (16S rRNA, dsrB and soxB) gradually decreased with increasing depth at S33 (Figure ) and were significantly negatively correlated with sediment depth. These results suggest that the sedimentary dynamic environment is likely an important factor in controlling the vertical distribution of the abundances of total bacterial 16S rRNA gene and functional genes (dsrB and soxB). Further analysis revealed similar vertical fluctuations in the abundance of the 16S rRNA, dsrB and soxB genes at S31, but the range of fluctuation obtained for the abundance of the 16S rRNA and soxB genes was smaller than that observed for the dsrB gene. In addition, the ratios of the functional genes (dsrB and soxB) to the total bacterial 16S rRNA gene were calculated: for the dsrB gene, the ratio ranged from 0.22 to 40.56% in S31 and from 0.06 to 4.2% in S33, and the ratios calculated for the soxB gene were less than 1% at both stations. These results indicate that the abundance of SRB might be more sensitive to the sedimentary dynamic environment than that of SOB and total bacteria.

The UPGMA clustering analysis results demonstrated different community structures of total bacteria, SRB, and SOB at the two sampling sites and among the different sediment depths, which might be due to the differences in the sediment depths and sedimentary dynamic environments between the two sampling stations. Overall, the sedimentary dynamic environment is most likely related to unmeasured variables, such as salinity, organic matter, dissolved oxygen, temperature, and H2S, which are widely accepted as important factors in controlling the community structures of total bacteria and sulfur-cycling bacteria (SRB and SOB) (Headd and Engel, 2013; Ye et al., 2016; Zhang et al., 2016). In addition, an RDA analysis revealed that sediment depth and DIN have a significant influence on the community structure of total bacteria, SRB, and SOB. In summary, our results suggest that environmental parameters, such as sediment depth and DIN, and the sedimentary dynamic environment, which showed differences between the two sampling stations, can significantly influence the community structure of total bacteria, SRB, and SOB.

Author Contributions

YZhe, TM, and ZY conceived and den class="Chemical">signed the experiments; YZha, XW, and HH performed the experiments and analyzed the data; YZha and YZhe wrote the paper.

Conflict of Interest Statement

The authors declare that the ren class="Gene">search was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Table 1

Vertical distribution of sediment bacterial community diversity and richness estimators based on the 16S rRNA gene.

Station	Depth (cm)	Read number	OTUs	Shannon	Chao 1	Good’s coverage (%)
S31	0–4	41,663	1,607	8.091	1770.117	98.9
	8–12	50,347	1,733	7.883	1975.441	98.5
	16–20	41,059	1,708	8.048	1951.907	98.7
S33	0–4	40,919	1,766	8.496	1818.074	99.2
	8–12	46,306	1,989	8.536	2237.209	98.6
	16–20	49,152	2,123	8.573	2325.423	98.5
	32–36	51,004	2,196	8.545	2640.667	98.2
	46–50	51,084	1,975	8.589	2297.087	98.6

Table 2A

Similarity-based OTUs and species richness and diversity estimates based on the dsrB gene.

Station	Depth (cm)	Read number	OTUs	Shannon	Chao 1	Good’s coverage (%)
S31	0–4	79,321	336	3.745	367.089	99.9
	8–12	78,861	240	3.782	266.464	99.9
	16–20	76,887	322	3.755	360.936	99.9
S33	0–4	67,284	336	2.708	343.483	99.9
	8–12	78,034	393	2.885	425.697	99.9
	16–20	76,646	456	3.247	485.94	99.9
	32–36	82,929	487	4.240	509.685	99.9
	46–50	75,100	413	4.27	426.448	99.9

Table 2B

Similarity-based OTUs and species richness and diversity estimates based on the soxB gene.

Station	Depth (cm)	Read number	OTUs	Shannon	Chao 1	Good’s coverage (%)
S31	0–4	12,100	180	5.060	185.571	99.8
	8–12	9,202	128	4.934	145.400	99.8
	16–20	9,070	103	4.057	115.720	99.8
S33	0–4	12,694	148	5.295	152.119	99.8
	8–12	12,817	141	4.845	142.747	99.9
	16–20	11,052	145	4.914	151.037	99.8
	32–36	18,786	113	5.086	110.563	99.9
	46–50	21,967	182	5.985	193.256	99.7

Table 3A

Rare and dominant OTUs based on the dsrB gene.

Station	Depth (cm)	Total number of OTUs	Dominant OTUs	Abundance of dominant OTUs (%)	Rare OTUs	Abundance of rare OTUs (%)
S31	0–4	336	13	83.54	183	0.78
	8–12	240	12	84.23	121	0.49
	16–20	322	12	81.41	180	0.70
S33	0–4	336	8	86.84	191	0.77
	8–12	393	9	86.70	257	1.03
	16–20	456	10	85.30	287	1.23
	32–36	487	11	77.16	220	0.97
	46–50	413	12	77.48	180	0.78

Table 3B

Rare and dominant OTUs based on the soxB gene.

Station	Depth (cm)	Total number of OTUs	Dominant OTUs	Abundance of dominant OTUs (%)	Rare OTUs	Abundance of rare OTUs (%)
S31	0–4	180	18	76.62	88	4.13
	8–12	128	19	78.45	47	1.76
	16–20	103	14	85.06	49	2.05
S33	0–4	148	23	76.96	61	2.53
	8–12	141	20	80.09	59	2.69
	16–20	145	21	76.78	59	2.51
	32–36	113	22	80.58	26	1.22
	46–50	182	30	67.58	50	1.72

Table 4A

SRB community composition based on taxonomic information of the dsrB and 16S rRNA genes.

Family level	OTUs		Total abundance (%)
	16S rRNA	dsrB	16S rRNA	dsrB
Peptococcaceae	6	6	0–0.33	0.45–6.60
Desulfobacteraceae	51	42	1.26–4.97	9.60–34.05
Desulfobulbaceae	38	4	0.27–9.37	0.11–5.67
Syntrophaceae	6	42	0.007–0.34	4.23–7.50
Nitrosomonadaceae	1	–	0.01–0.62	–
Desulfarculaceae	26	–	0.32–5.82	–
Desulfovibrionaceae	1	–	0–0.01	–
Syntrophobacteraceae	52	–	0.13–0.77	–
Nitrospiraceae	12	–	0.10–1.16	–
Thermodesulfovibrionaceae	42	–	1.72–8.59	–

Table 4B

SOB community composition based on taxonomic information of the soxB and 16S rRNA genes.

Family level	OTUs		Total abundance (%)
	16S rRNA	soxB	16S rRNA	soxB
Spirochaetaceae	30	8	0.03–0.69	0.85–5.04
Chlorobiaceae	–	2	–	0.13–13.08
Helicobacteraceae	5	–	0.02–4.57	–
Rhodospirillaceae	42	5	0.24–1.43	0.12–3.75
Rhodobacteraceae	27	20	0.27–3.32	1.33–14.74
Bradyrhizobiaceae	1	35	0–0.01	15.79–36.26
Hyphomicrobiaceae	16	2	1.38–2.67	0.01–3.51
Rhizobiaceae	1	–	0–0.01	–
Neisseriaceae	1	–	0–0.03	–
Comamonadaceae	7	3	0.02–0.15	0.05–2.88
Ectothiorhodospiraceae	10	2	0.07–0.56	0.09–6.38
Thiotrichaceae	7	2	0–0.07	0.06–1.60
Piscirickettsiaceae	44	4	0.52–13.67	0–3.38
Burkholderiaceae	–	6	–	0.13–4.82
Oxalobacteraceae	–	1	–	0–4.36
Hydrogenophilaceae	–	1	–	0–0.44
Chromatiaceae	–	18	–	0.25–7.39
Halothiobacillaceae	–	2	–	0–1.42