Thirty-four NHFS (N1, N2, N7-N29, MS1, MS2, and M1-M7) and five HFS samples (N3-N6 from Wocan and M8 from Longqi) were collected along the CR (N1-N29) and SWIR (MS1, MS2, and M1-M8) ( Figure 1 ) during the DY115-21, DY115-22, DY125-26, and DY-28I cruises of the DaYang YiHao and Haiyang 20 research vessels (RVs) using TV-guided grabs or multi-core samplers (Φ10 × 60 cm). The samples were transferred into sterile 200-ml plastic boxes and frozen at -80°C until further processing in the laboratory. The longitude, latitude, water depth, and temperature of the sampling sites were recorded by the RVs.

The physicochemical analysis of the samples was conducted at the Analysis Center of Tsinghua University using the methods as previously described (Chen et al., 2016). Detailed information of the samples was provided in Table S1 . Significant differences of environmental factors between the HFS and NHFS samples were determined using the Wilcoxon test. The map of sampling sites was drawn using the Generic Mapping Tools (GMT) (Wessel and Smith, 1998) (http://gmt.soest.hawaii.edu/) and plotted using the Mercator projection and NUVEL-1A plate boundaries (http://jules.unavco.org/GMT/).

The total DNA was extracted from 5 g (wet weight) of each of the 39 sediment samples which was carefully homogenized before subsampling. The procedure was as previously described (Chen et al., 2016). Amplifications were conducted using the primer set F515 (5’-GTGCCAGCMGCCGCGGTAA-3’) and R806 (5’-GGACTACHVGGGTWTCTAAT-3’) that targets V4 region of the prokaryotic 16S rRNA genes. The primers were labeled with a unique 6-nt barcode for each sample. The PCR of each sample was performed in triplicate with a 25 μL reaction mixture containing 10 ng purified DNA template under following thermal cycles: initial denaturation at 98°C for 1 min; 35 cycles at 98°C for 10 sec, 50°C for 30 sec, and 72°C for 30 sec; and a final extension at 72°C for 5 min. The triplicate PCR products were mixed together and purified using the Takara PCR fragment purification kit (Takara, Japan), followed by the ligation of ‘A’ sequencing adaptor. The purified PCR amplicons were combined in equimolar amounts and submitted to Novogene Co. Beijing for Illumina 250 bp paired-end sequencing on a MiSeq machine.

Fragments (~350 bp) of bacterial dsrB gene related to the sulfur cycle was amplified using primers DSRp2060F (5’-CAACATCGTYCAYACCCAGGG-3’) and DSR4R (5’-GTGTAGCAGTTACCGCA-3’) (Wagner et al., 1998; Geets et al., 2006). The PCR was performed in triplicate with a 50 μL reaction mixture containing 50 ng purified DNA template under following thermal cycles: initial denaturation at 94°C for 2 min; 15 cycles at 94°C for 10 sec, 62°C for 30 sec, and 72°C for 20 sec; 20 cycles at 94°C for 10 sec, 48°C for 15 sec, and 72°C for 20 sec; and a final extension at 72°C for 2 min. The PCR products were purified using the Wizard SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA), and then submitted to Beijing Genomics Institute (BGI) Tech., where the library was constructed and sequenced using Illumina Miseq PE 300.

Raw fastq sequences were first demultiplexed according to the unique barcode by the software QIIME v 1.9.1 (Quantitative Insights Into Microbial Ecology) (Caporaso et al., 2010b), and then pairs of reads were merged with FLASH (Fast Length Adjustment of Short Reads) (Magoč and Salzberg, 2011) under default settings. The merged reads were quality filtered by QIIME. Sequences that contained ambiguous bases (N) were discarded, and sequences were truncated at any sites containing more than three consecutive bases with a Phred quality score (Q) < 20; only those with > 75% (as percentage of total read length) consecutive high-quality base calls were included in further analyses (Bokulich et al., 2013). Chimeric sequences were detected and excluded using the QIIME implementation of the UCHIME algorithm (Edgar et al., 2011) based on the Greengenes database version 13.8 (DeSantis et al., 2006).

Statistics analyses of the high-quality 16S rRNA gene sequences were performed with QIIME and USEARCH as described previously (Chen et al., 2016). Briefly, the 16S rRNA gene sequences were assigned to operational taxonomic units (OTUs) at a similarity cutoff of 97% by UPARSE (Edgar, 2013). After discarding all singletons, representative sequences of the OTUs were aligned using PyNAST (Caporaso et al., 2010a). The phylogenetic tree of the representative sequences was constructed using FastTree (Price et al., 2010). The RDP classifier (Cole et al., 2014) was used to annotate taxa based on the SILVA database version 132 (Yilmaz et al., 2014) with a 50% bootstrap cutoff. To avoid sampling errors, subsamples of the minimal number of sequences (9,231 sequences) per sample were randomly selected to compare the alpha and beta diversity between samples. The estimators of Shannon diversity, Chao1 richness, and Faith’s phylogenetic diversity (PD) were calculated in QIIME. The Wilcoxon test was used to measure the difference of environmental factors between the HFS and NHFS samples.

Bray-Curtis distance and Euclidean distance matrices were used to measure the dissimilarities of prokaryotic communities and geochemical characteristics, respectively, and the non-metric multidimensional scaling (NMDS) analysis was used to visualize the dissimilarities between the samples. Nonparametric analysis of similarities (ANOSIM) was performed on Bray-Curtis community dissimilarities. Standardization was performed to make all geographic and physicochemical parameters comparable between the different units of different variables using the scale function in R v4.0.2 (https://www.r-project.org/). The “Geography” library was used to calculate the paired geographic distances between samples based on latitude and longitude coordinates in R. A Spearman’s rank correlation between the communities and geochemical factors was calculated. Mantel test, Redundancy analysis (RDA), and Bioenv analysis were conducted to assess the relationships between prokaryotic communities and environmental variables with 999 permutations. These analyses were carried out with the vegan package v2.5-7 in R (Oksanen et al., 2020). The differentially abundant taxa were determined by the linear discriminant analysis (LDA) effect size (LEfSe) online software (http://huttenhower.sph.harvard.edu/lefse/) (Segata et al., 2011). The threshold on the logarithmic LDA score for distinguishing features was set to 3.0. The software PICRUSt2 (Douglas et al., 2020) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database were used to infer functional potentials of the microbial communities based on the 16S rRNA gene amplicon sequencing profile. STAMP v2.0.0 (Parks et al., 2014) was used to determine significant putative KEGG pathways with two-side Welch’s t-test (P value < 0.05) and the Benjamini-Hochberg FDR (Benjamini and Hochberg, 1995) correction for multiple tests.

The sequence processing was the same as above. The USEARCH software (Edgar, 2010) was used to divide the dsrB sequences into different OTUs based on 90% sequence similarity, removing the chimera by the denovo method. All OTUs were annotated using the NCBI blastx method (Altschul et al., 1997) and Non-redundant protein sequence (NR) database (Pruitt et al., 2007). After discarding the sequences with E-value > e^-10, coverage < 50%, or bit-score < 50, and filtering out the false positive sequences annotated as non-dsrB genes, the high-quality sequences were obtained for subsequent analysis. Alpha diversity of the dsrB sequences was analyzed as above. The sequences were compared to those of the DsrAB database constructed by Müller (Müller et al., 2015) through local blast analysis to retrieve the most similar reference sequences. The phylogenetic tree of the representative dsrB sequences and related reference sequences was constructed using MEGA7 (Kumar et al., 2016).

The metadata of the 16S rRNA gene and dsrB amplicon sequences from this study have been deposited in the NCBI Sequence Read Archive (SRA) under the accession numbers SRR17439176-SRR17439214 and SRR17440762-SRR17440769, respectively, and in the China National Microbiology Data Center (NMDC) under the accession numbers NMDC40016400-NMDC40016438 and NMDC40016793-NMDC40016800, respectively.

Three geographic and 12 physicochemical parameters of the samples were obtained and are summarized in Table S1 . Five samples were from HFSs at water depths of 2800-3104 m and with pairwise geographic distances ranging from 0.1 to 5049 km; 34 samples were from NHFSs at water depths of 1608-4213 m, and the pairwise geographic distances between the sampling sites ranged from 1 to 5582 km ( Figure 1 ). NMDS analysis based on the 12 physicochemical factors and water depth was used to evaluate the dissimilarities between the samples. As shown in Supplementary Figure S2 , HFSs were obviously different from NHFSs while most SWIR samples were similar to CR samples, indicating that the environmental variables were significantly affected by habitat rather than by geographic distance. Compared to NHFSs, HFSs contained significantly lower content of Ca (Wilcoxon test, P < 0.0001) and lower pH (P < 0.001) while higher contents of TP (P < 0.001), Fe, S and TOC (all P < 0.05) ( Figure 2 ).

After quality control and chimera removal, a total of 1,534,540 16S rRNA gene sequences were obtained from the 39 samples, including 1,141,817 bacterial and 392,723 archaeal sequences. These sequences were binned into 14,978 non-singleton OTUs (1,663-6,260 OTUs per sample) at a 97% sequence similarity cutoff, with a coverage of 84.75-92.31% across the samples ( Table S2 ). The OTUs belonged to 74 prokaryotic phyla ( Supplementary Figure S3 ), 235 classes and 2,059 genera (892 classified genera and 1,167 unclassified genera). Proteobacteria (relative abundance 36.98%), Thaumarchaeota (23.31%), Planctomycetes (8.53%), Actinobacteria (5.44%), Acidobacteria (5.33%), Chloroflexi, Bacteroidetes, Gemmatimonadetes, and Nitrospirae were the dominant phyla (average relative abundance >1% in all samples) in the communities. The dominant classes (relative abundance >1%) were Nitrososphaeria, Gammaproteobacteria, Alphaproteobacteria, Deltaproteobacteria, Acidimicrobiia, Phycisphaerae, Planctomycetacia, Bacteroidia, Dehalococcoidia, Gemmatimonadetes, Nitrospira, and Anaerolineae, constituting 78.75% of the total sequences. At the genus level, Nitrosopumilus of the phylum Thaumarchaeota was the most abundant (7.31%), followed by unclassified_Nitrosopumilaceae (5.67%), Woeseia (3.95%), Escherichia-Shigella (2.85%), Urania-1B-19 marine sediment group (2.68%), unclassified_Actinomarinales (2.17%), unclassified_Kiloniellaceae (1.69%), Halomonas (1.63%), Sphingomonas (1.61%), Nitrospira (1.56%), uncultured_Gemmatimonadaceae (1.55%), and Pir4 lineage (1.31%).

There was no significant difference in alpha diversity of the community between the HFS and NHFS habitats or between the CR and SWIR locations (P > 0.05; Supplementary Figure S4 ). However, NMDS analysis at 97% OTU similarity of 39 subsamples of even reads showed a biogeographic distribution of the microbial community between HFSs and NHFSs but not between the CR and SWIR ( Figure 3A ). This result accords with the distribution of environmental variables ( Supplementary Figure S2 ).

Most of the above nine dominant phyla showed remarkably different relative abundances between HFSs and NHFSs ( Figure 3B ). Besides, Firmicutes, Nanoarchaeaeota, and Patescibacteria were also dominant in the HFS samples ( Figure 3B ). Using LEfSe at the phylum level, 16 of the 74 phyla were determined to be differentially abundant between HFSs and NHFSs (LDA > 3.0, P < 0.05; Figure 4 ). Nine phyla were enriched in HFSs, of which Proteobacteria, Aerophobetes, Bacteroidete, and Deferribacteres were the most enriched (LDA > 4.0, P < 0.05; Figure 4 ); seven phyla were enriched in NHFSs, with Thaumarchaeota the most enriched (LDA > 4.0, P < 0.05; Figure 4 ). It should not be neglected that Epsilonbacteraeota and Hydrothermarchaeota, previously reported to be enriched in HFSs (Kato et al., 2010; Zhou et al., 2020), were also enriched in our HFS samples (LDA > 3.0, P < 0.05; Figure 4 ). Epsilonproteobacteria are important chemolithotrophic primary producers in deep-sea hydrothermal vent systems (Kato et al., 2010; Waite et al., 2017), and Hydrothermarchaeota are widely distributed in HFSs (Zhou et al., 2020). Among the 2,059 genera, 54 (16 classified and 38 unclassified genera, accounting for 54.36% of all reads) were detected to be differentially abundant between HFSs and NHFSs by LEfSe (LDA > 3.0, P < 0.05; Supplementary Figure S5 ). HFSs and NHFSs were enriched for 28 and 26 genera, respectively. Achromobacter and Sphingomonas of Proteobacteria, uncultured_Bacteroidetes and Ambiguous_Actinobacteria were the most enriched in HFSs; while Nitrosopumilus and three uncultured genera of Thaumarchaeota and the uncultured Urania-1B-19 group of Planctomycetes were the most enriched in NHFSs (LDA > 4.0, P < 0.05, Supplementary Figure S5 ). Members of the genera Achromobacter and Sphingomonas have previously been found in deep-sea sediments (Zhang et al., 2014; Yang et al., 2020b; Kikukawa et al., 2021), and this is the first study to demonstrate that they are enriched in HFSs. In the previous studies of microbial diversity in deep-sea hydrothermal vents, Nitrosopumilus was found to adapt to inactive vents or cooler ecological niches of the vents (Zhang et al., 2016; Cerqueira et al., 2017). Our study further suggests that Nitrosopumilus might prefer NHFSs. As Nitrosopumilus consists of marine ammonia-oxidizing archaea, this genus may play a significant role in the carbon and nitrogen cycling in non-hydrothermal environments of the global ocean (Walker et al., 2010; Qin et al., 2017).

Although the microbial community compositions were obviously different between HFSs and NHFSs, 190 OTUs were commonly present in all samples. These common OTUs accounted for only 1.27% of the total OTUs but 50.81% of the total reads, with the relative abundance of each OTU ranging from 0.02 to 4.93%. These OTUs belonged to 13 phyla ( Supplementary Figure S6A ), most of which have been reported to be common phyla in marine. Only 11 of these OTUs had a relative abundance of >1% each, and seven of them belonged to Nitrosopumilaceae of Thaumarchaeota. The phyla Proteobacteria (13.56%) and Thaumarchaeota (21.02%) were the most dominant in the common OTUs of HFSs and NHFSs, respectively ( Supplementary Figure S6B ). At the genus level, these common OTUs belonged to 106 genera, and 78 of them were unclassified, indicating that most of them represent novel taxa that have been little studied.

Our data on the SWIR samples contained all the phyla that were detected in a previous clone library study on semi-consolidated carbonate sediments of the SWIR (Li et al., 2014) ( Supplementary Figure S1 ). The relative abundances of most of these phyla were largely consistent between the studies, but Planctomycetes, Chloroflexi, Nitrospirae, and Thaumarchaeota were more abundant while Acidobacteria and Verrucomicrobia were less abundant in our study. In addition, one dominant phylum detected here, Gemmatimonadetes, was absent in the previous study. The dominant prokaryotic taxa showed in HFSs of our study were contrast to those reported in a recent study on the newly discovered serpentinite-hosted Old City hydrothermal field of the SWIR (Lecoeuvre et al., 2021). The Old City hydrothermal chimneys had low relative abundances of Actinobacteria, Bacteroidetes, Proteobacteria, and Nanoarchaeaeota while more abundant Schekmanbacteria and Thaumarchaeota. The varied microbial communities in different samples of the SWIR may be attributed to different local geochemistry factors.

The concordant NMDS biplots of prokaryotic communities and physicochemical factors suggested a linkage between the community structure and environmental factors ( Figure 3A and Supplementary Figure S2 ). This relationship was supported by the significant Spearman correlation between the community compositions and the environmental attributes of all 39 samples ( Figure 5A ). RDA analysis revealed that all of the environmental factors explained as much as 44.85% of the overall variance of the communities (P = 0.001), with Ca, S, TP, and TN showing significant effects ( Figure 5B ; Table S3 ). Bioenv analysis indicated that water depth, Ca, S, nitrate nitrogen, Fe, and TN were the main environmental variables affecting the beta diversity (dissimilarity) among communities, explaining as much as 36.04% of the overall variance. As the content of Ca, S, TP, and Fe differed significantly between HFSs and NHFSs ( Figure 2 ), these results sustain that the environmental differences between the two habitats shape the prokaryotic community structure. The results also align with our previous finding that Ca and TP, as essential nutrients, significantly correlate to prokaryotic communities along the SWIR (Chen et al., 2016). Moreover, S and Fe have been reported to be enriched in deep-sea hydrothermal vents (Jannasch and Mottl, 1985) and play an important role in fueling the microbes present in the vent habitats on the SWIR (Cao et al., 2014).

To quantify the contribution of spatial distance to prokaryotic community variability, Spearman correlations between Bray-Curtis community dissimilarities and geographic distances were calculated. There was no significant correlation across all samples. However, when focusing on the NHFS samples only, a weak significant correlation was observed (r = 0.195, P = 0.022; Supplementary Figure S7 ). The NHFS samples had relatively similar environmental variables but geographic distance up to 5582 km, it was therefore expected that dispersal limitation would become important at such a spatial scale. This is consistent with the commonly observed distance–decay pattern (Martiny et al., 2011). For the whole dataset containing both the HFS and NHFS samples, the remarkable difference in environmental variables probably disrupts the overall distance–decay pattern, supporting the Baas-Becking hypothesis that ‘the environment selects’ (Martiny et al., 2006). Together, the prokaryotic biogeographic pattern along the IOR is readily shaped by environmental selection between HFSs and NHFSs while also affected by large spatial distance (>5000 km).

To explore the functional capabilities of the prokaryotic communities in the sediments of the IOR, we used PICRUSt2 (Douglas et al., 2020) to predict the metabolic functions from the 16S rRNA gene sequences. A total of 437 metabolic pathways (level 3) belonging to 49 level 2 and seven level 1 KEGG pathway groups were predicted in all the samples. As shown in Supplementary Figure S8 , nine level 2 KEGG pathway groups were significantly different in abundance between the HFS and NHF communities (P < 0.05). The HFS communities were mainly enriched for pathways related to lipid metabolism, unclassified: signaling and cellular processes, and cell growth and death ( Supplementary Figure S8 ), reflecting the ability of microbes to adapt to dynamic changes in hydrothermal habitats. On the other hand, pathways of energy metabolism, metabolism of cofactors and vitamins, and transcription were mainly enriched in the NHFS communities. The microbes in NHFS may be committed to obtaining more energy to maintain physiological metabolic activities.

Energy metabolism pathways related to the cycles of carbon, nitrogen, and sulfur elements showed some differences between the HFS and NHFS samples ( Figure 6 ). Compared to HFSs, NHFSs contained more abundant genes involved in carbon fixation pathways (P = 0.021; Figure 6A , Supplementary Figure S9 ). The genes required for the 3-hydroxypropionate/4-hydroxybutyrate (3-HP/4-HB) pathway for CO₂ fixation, i.e., accA encoding the acetyl-coenzyme A (CoA) carboxylase and abfD encoding the 4-HB-CoA dehydratase, were predicted. Although accA was abundant in HFSs, abfD was lacked in the samples; whereas both genes were abundant in NHFSs ( Figure 6B ). The abundant Thaumarchaeota and some Crenarchaeota in NHFSs likely use the 3-HP/4-HB pathway for autotrophic carbon fixation (Pester et al., 2011). Meanwhile, the three key genes of the reductive tricarboxylic acid (rTCA) cycle, namely aclB, oorA, and por/nifJ, were also predicted with a higher abundance in NHFSs. Microbes utilizing the rTCA cycle to fix CO₂ are ubiquitous in microaerophilic to anaerobic and sulfuric-rich environments, such as deep-sea hydrothermal vents (Campbell and Cary, 2004; Nakagawa and Takai, 2008). As for methane metabolism pathways, the genes encoding copper-dependent particulate methane monooxygenase (pMMO) for methane oxidation (pmoA, pmoB, and pmoC) were enriched in NHFSs, probably contributed by methane-oxidizing bacteria of Alphaproteobacteria, Gammaproteobacteria, and NC10 (Kirchman, 2012; Offre et al., 2013). Moreover, methanotrophs seemed to overwhelm methanogens in both HFSs and NHFSs, as the average abundance of methane-oxidizing genes such as pmo were more than 10 times that of methanogenic genes such as mcr. This may be attribute to the fact that methane is abundant in deep sea sediments (Baker et al., 1995).

A number of gene families related to nitrogen metabolism were predicted ( Figure 6 ), forming a complete nitrogen cycle. In general, nitrogen metabolism pathways were enriched in NHFSs (P = 0.048) ( Figure 6A ). The gene families involved in N₂-fixation, including nifD, nifH, nifK, and nifW, showed little difference in relative abundance between HFSs and NHFSs ( Figure 6B ). It was noticeable that the genes for ammonia oxidation (amoA, amoB, and amoC) were much more abundant in NHFSs than in HFSs ( Figure 6B ). Congruently, prokaryotes related to ammonia oxidation, especially members of Thaumarchaeota, had a relative abundance in NHFSs more than 10 times that in HFSs ( Figure 4 and Supplementary Figure S5 ). These results support that N cycle is more important in NHFSs than in HFSs, because nitrification, containing ammonia oxidation and nitrite oxidation, is the important process in N cycle and ammonia oxidation is the rate-limiting step of nitrification. Furthermore, gene families for each process of denitrification were predicted, including napAB, nirKS, norBC, and nosZ, where nirK (encoding Cu-containing nitrite reductase, Cu-NIR) was remarkably more abundant in NHFSs than in HFSs. Meanwhile, some assimilatory and dissimilatory nitrate reduction genes (nasA, nirBD, and nrfA) were enriched in HFSs.

Collectively, the gene families involved in carbon fixation and ammonia oxidation pathways were more abundant in NHFSs than in HFSs. Both of the two pathways seem to be related to Thaumarchaeota, which was significantly enriched in NHFs according to the result of LEfSe ( Figure 4 ). It has been illustrated that some members of Thaumarchaeota support their chemolithoautotrophic lifestyle in the process of ammonia oxidation through the 3-HP/4-HB pathway (Pester et al., 2011). Given the high relative abundances of both the genes and taxa participating in the N and C cycles in NHFSs, the N cycle might be centralized in an energy-generating manner to support the microbial community in NHFSs.

On top of that, diverse genes involved in the S cycle were also predicted ( Figure 6B ). The assimilatory sulfate reduction pathway was abundant compared to the other pathways in both HFSs and NHFSs. Assimilatory sulfate reduction is the predominant way for microbes to integrate sulfur from environment into cells, and plays a vital role in microbial sulfur metabolism. On one hand, the sulfate assimilation process is the key source of cysteine in microorganisms; on the other hand, this pathway also involves in the detoxification of heavy metals in the organisms (Salazar et al., 2013; Wheaton et al., 2016), which may help deep-sea microbes to withstand high concentrations of metal ions. Sulfate represents one of the most common electron acceptors in the ocean. Genes related to the dissimilatory sulfate reduction pathway were mainly found from Alphaproteobacteria, Betaproteobacteria, Deltaproteobacteria, and Gammaproteobacteria in the samples and assigned to the apr, dsr, and qmo families. These gene families had a slightly higher relative abundance in HFS than in NHFS ( Figure 6B ). The dsr genes might be involved in both dissimilatory sulfur reduction and oxidation. Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria were also the major classes related to sulfur oxidation, contributing gene families such as soxABCXYZ, fccAB, sqr and sseA. To sum up, the S cycle-related pathways were relatively complete and enriched in both the HFS and NHFS samples with no significant difference between the two habitats.

Although the overall relative abundance of sulfur metabolism pathways was not significantly different between HFSs and NHFSs ( Figure 6 ), prokaryotes with sulfate-reducing ability exist widely in deep-sea hydrothermal vents (Frank et al., 2013; Meier et al., 2019). To further explore the difference of sulfate-reducing prokaryotes between HFSs and NHFSs, we used four HFSs (N3, N4, N5, and N6) and four NHFSs (N8, N12, N15, and N18) from CR as representative samples to investigate the diversity and distribution of dsrB genes, the well-known molecular marker of sulfate reducers. NMDS analyses of the physicochemical factors and 16S rRNA gene-based microbial communities both showed that the eight samples also clustered according to habitat ( Supplementary Figure S10 ). A total of 427,499 high-quality dsrB sequences were obtained, which were divided into 129 OTUs according to 90% nucleotide similarity. The two habitats shared 25 dsrB OTUs, and the HFS samples showed a significant higher alpha diversity of dsrB than the NHFS samples (P < 0.05; Table S4 ). A comparison to the DsrAB reference database (Müller et al., 2015) and subsequent phylogenetic analysis ( Figure 7 ) showed that the OTUs were assigned into 13 class-level taxa, mainly uncultured and unclassified clades. The vast majority of dsrB OTUs (92/129) were clustered with sequences from marine sources. Only a small proportion of the OTUs were clustered with sequences of known sulfate-reducing strains, most of which belong to Deltaproteobacteria, dominated by Desulfobulbaceae and followed by Desulfobacteraceae and Syntrophobacteraceae. Very few OTUs were affiliated to Clostridia of the phylum Firmicutes. Deltaproteobacteria and Firmicutes are already known as common sulfate-reducing bacteria (Kaneko et al., 2007; Robador et al., 2014). In addition, dsrB sequences related to Thermodesulfobacteria (OTU35) were detected only in the HFS samples, in line with the thermophilic characteristic of the taxon. A previous study found that Thermodesulfobacteria existed in high abundance as sulfate-reducing microorganisms in deep-sea hydrothermal chimneys in the Okinawa trough and Izu-Bonin arc (Nakagawa et al., 2004). Forty-two OTUs (61.7% of the total sequence) had a similarity of < 90% with the known sequences in the reference (Müller et al., 2015) and NR database, and thus could be considered as putative novel dsrB genes. Among them, OTU86 was the most dominant OTU (26.3% of the total sequence), existing in all samples and showing a low similarity (< 77%) to the known dsrB sequences. This OTU may represent a potential novel bacterial phylogroup widely distributed in CR sediments.

In short, dsrB sequencing analysis further confirms the existence of diverse and novel sulfate/sulfite-reducing microbes in the IOR sediments. The HFS samples contain higher diversity of sulfate reduction genes and taxa than the NHFS samples, inferring that the S cycle is more important for the HFS microbial community.