Distinct Bottom-Water Bacterial Communities at Methane Seeps With Various Seepage Intensities in Haima, South China Sea

Methane seeps are chemosynthetic ecosystems in the deep-sea environment. Microbial community structures have been extensively studied in the seepage-affected sediments and investigation in the water column above the seeping sites is still lacking. In this study, prokaryotic communities in the bottom water about 50 cm from the seabed at methane seeps with various seepage intensities in Haima, South China Sea were comparatively studied by using 16S ribosomal RNA gene sequencing. These sites were assigned based on their distinct methane content levels and seafloor landscapes as the non-seepage (NS) site, low-intensity seepage (LIS) site, and high-intensity seepage (HIS) site. The abundances of the dominant phyla Proteobacteria, Bacteroidetes, and Actinobacteria differed significantly between NS and the two seepage sites (p < 0.05). Alpha diversity differed among the three sites with the HIS site showing the lowest community diversity. Principal component analysis revealed highly divergent bacterial community structures at three sites. Many environmental variables including temperature, alkalinity, pH, methane, dissolved organic carbon (DOC), and inorganic nutrients were measured. Redundancy analysis indicated that methane content is the key environmental factor driving bacterial community variation (p = 0.001). Linear discriminant analysis effect size analysis identified various differentially enriched genera at the LIS and HIS sites. Phylogenetic analysis revealed close phylogenetic relationship among the operational taxonomic units of these genera with known oil-degrading species, indicating oil seepage may occur at the Haima cold seeps. Co-occurrence networks indicated that the strength of microbial interactions was weakest at the HIS site. This study represents a comprehensive comparison of microbial profiles in the water column of cold seeps in the SCS, revealing that the seepage intensity has a strong impact on bacterial community dynamics.

Methane seeps are chemosynthetic ecosystems in the deep-sea environment. Microbial community structures have been extensively studied in the seepage-affected sediments and investigation in the water column above the seeping sites is still lacking. In this study, prokaryotic communities in the bottom water about 50 cm from the seabed at methane seeps with various seepage intensities in Haima, South China Sea were comparatively studied by using 16S ribosomal RNA gene sequencing. These sites were assigned based on their distinct methane content levels and seafloor landscapes as the non-seepage (NS) site, low-intensity seepage (LIS) site, and high-intensity seepage (HIS) site. The abundances of the dominant phyla Proteobacteria, Bacteroidetes, and Actinobacteria differed significantly between NS and the two seepage sites (p < 0.05). Alpha diversity differed among the three sites with the HIS site showing the lowest community diversity. Principal component analysis revealed highly divergent bacterial community structures at three sites. Many environmental variables including temperature, alkalinity, pH, methane, dissolved organic carbon (DOC), and inorganic nutrients were measured. Redundancy analysis indicated that methane content is the key environmental factor driving bacterial community variation (p = 0.001). Linear discriminant analysis effect size analysis identified various differentially enriched genera at the LIS and HIS sites. Phylogenetic analysis revealed close phylogenetic relationship among the operational taxonomic units of these genera with known oil-degrading species, indicating oil seepage may occur at the Haima cold seeps. Co-occurrence INTRODUCTION Cold seeps are highly productive ecosystems dominated by profuse chemosynthetic microbes (Boetius and Wenzhofer, 2013;Joye, 2020). Apart from thermal maturation of organic matter at depth (Galimov, 1988), methane (CH 4 ) is also produced through degradation of organic matter in the subsurface sediment by methanogens; sulfate-reducing bacteria and anaerobic methanotrophic (ANME) archaea consume CH 4 in the sediments that is often anoxic and aerobic methanotrophic bacteria in the oxic water column oxidize CH 4 seeped from sediment (Niu et al., 2018). These processes facilitate the formation of carbonates and result in high concentrations of hydrogen sulfide and oxygen depletion . In general, these geochemical alterations may help to reconstruct the microbial communities in the waters above cold seeps. Chemosynthetic organisms in cold seeps utilize CH 4 or hydrogen sulfide to sustain various faunal assemblages (Pasulka et al., 2016) such as tubeworms (Zhao et al., 2020), mussels (Sun et al., 2017;Xu et al., 2019), and bivalves.
Because abundant cold seeps are distributed in the South China Sea (SCS), it is received increasing research interest. To date, more than 30 cold seeps have been discovered on the northern continental slope of the SCS , but only two cold seeps are still active and identified as harboring faunal assemblages . One is site F and another is Haima seeps, which reported for the first time in 2015, located at a depth of approximately 1,360-1,400 m in the southern part of the Qiongdongnan Basin .
Although microbial community structures have been investigated at cold seeps in the SCS, nearly all the reports have focused on seepage-affected sediments (Lai et al., 2007;Zhang et al., 2012;Niu et al., 2017;Wu et al., 2018;Cui et al., 2019Cui et al., , 2020Zhuang et al., 2019;Li et al., 2020). Recently, Zhang et al. (2020) investigate prokaryotes along the depth profile of the water column. These authors found heterotrophs to be the major component of the bacterial community with the most abundant functions being enriched in carbohydrate metabolism in samples from deeper water layers of the Haima seep, indicating the presence of ample organic matter. In addition, Guan et al. (2018) detected unresolved complex mixtures (UCMs) in the hydrocarbon fractions of authigenic carbonates retrieved from Haima seeps, implicating the presence of heavily biodegraded organic matter or, more likely, oil. Indeed, UCMs have frequently been observed in authigenic seep carbonates from sites where oil seepage occurs (Birgel et al., 2011;Feng et al., 2014). Apart from UCMs, lipid biomarker patterns in the authigenic carbonates also indicate oil seepage in the Haima cold seeps (Guan et al., 2021). In addition, mud diapirs and faults are well developed in the Haima cold seeps, potentially serving as pathways for the upward migration of gas fluids from deep reservoirs to the seafloor (Xie et al., 2006;Huang et al., 2016). Therefore, an oil reservoir may exist in the Haima cold seeps and microbial communities in the water column may have the ability to degrade oil.
Although the microbial community in the water column of seep sites in the SCS has been studied, a comparison of sites with various seepage activities is still lacking. To elucidate how seepage intensity shapes the prokaryotic community, we compared the prokaryotic communities of seep sites in Haima with distinct seafloor landscapes and CH 4 contents. Bottomwater samples were collected from sites of the non-seepage (NS), low-intensity seepage (LIS), and high-intensity seepage (HIS). Bacterial and archaeal communities were studied by using 16S ribosomal RNA (rRNA) gene amplicon sequencing and the concentrations of geochemical factors in the seep environment were quantified. The main aims were to better understand: (i) whether intergroup variance exists among sampling sites with different seepage intensities, (ii) the correlations between environmental factors and the prokaryotic community structure, and (iii) the dominant/signature microbes in active cold seeps, especially with regard to whether oil degraders inhabit such ecosystems.

Collection of Bottom-Water Samples and Geochemical Analysis
Sampling was conducted during an expedition in May 2019 at the Haima cold seeps located in the SCS, China. Bottom-water samples were collected by using 10-L Niskin bottles fixed on Remote Operated Vehicle (ROV) Haima from depths of 1,360-1,398 m ( Figure 1A). Three sample sites were selected based on distinct seafloor landscapes: a NS site outside active seepage (Figure 1B), the LIS site with dead clams and live tubeworms (Figure 1C), and the HIS site with live mussels and continuous bubbling of gas ( Figure 1D). Six samples from six stations were collected from the NS and HIS site and four samples were collected from the LIS site. The distance between three sites was 1-3.7 km and the stations at each site were approximately 100 m apart between each other.
In situ CH 4 concentrations were measured at each site by using the HydroCTM Methane Sensor (HYDROC GmbH, Flensburg, Germany, United Kingdom) mounted on the ROV Haima. In addition, the concentrations of inorganic nutrients, including nitrate, nitrite, ammonia, and phosphate, and the pH were quantified in reference to those described in "Specification of Oceanographic Investigation in China" (GB/T12763.4-2007, General Administration of Quality Supervision, Inspection, and Quarantine of the People's Republic of China, 2008). 1 The alkalinity was determined on board by direct titration with ∼0.006 M hydrochloric acid (HCl) by using a pH meter. A total of 4 L of collected liquid were immediately filtered through a 0.22µm polycarbonate membrane (47-mm diameter, Millipore, MA, United States) to collect microbes. All the filtrates were stored at -80 • C until further analysis.

Deoxyribonucleic Acid Extraction, PCR Amplification, and Sequencing
Total genomic DNA of microbes was extracted from the 0.22µm polycarbonate filters by using the ALFA-SEQ Advanced Water DNA Kit (Magigene, Guangdong, China) following the instructions of the manufacturer. The extracted DNA was 1 http://www.gb688.cn/bzgk/gb/newGbInfo?hcno= D7C5F44155DBE2C0BC40EA764A6BBF4A quantified with a NanoDrop 2000 Spectrophotometer (Thermo Scientific, Thermo Fisher Scientific Corporation, Waltham, MA) and the quality was checked via gel electrophoresis. The bacterial 16S rRNA gene V3-V4 region was amplified by using the specific forward primer 314F 5 -CCTAYGGGRBGCASCAG-3 and reverse primer 806R5 -GGACTACNNGGGTATCTAAT-3 (Yu et al., 2005). The Archaea 16S rRNA gene V4 region was amplified with the forward primer Uni519F 5 -CAGYMGCCRCGGKAAHACC-3 and reverse primer Arch806R 5 -GGACTACNSGGGTMTCTAAT-3 (Zhang et al., 2015). Sample-specific 12-bp barcodes were assigned to perform multiplex sequencing. PCR was carried out by using the TaKaRa Premix Taq R version 2.0 (TaKaRa Biotechnology Corporation, Dalian, China) containing 25 µL premix Taq (2X), 1 µL (10 µM) each forward and reverse primer, 50 ng DNA template, and deionized distilled water (ddH 2 O) to a total of 50 µL. The "contaminant" control was performed with sterile ddH 2 O as the DNA template. PCR amplification was conducted as follows: 94 • C for 5 min followed by 30 cycles consisting of denaturation at 94 • C for 30 s, annealing at 52 • C for 30 s, and extension at 72 • C for 30 s and a final extension of 10 min at 72 • C. PCR amplicons from PCR triplicates were mixed and examined by 1% agarose gel electrophoresis. The amplicons of each sample were pooled in equal amounts after concentration measurement with the GeneTools Analysis Software (version 4.03.05.0, SynGene, Karnataka, India) and purified with the E.Z.N.A. R Gel Extraction Kit (Omega, United States). Library construction was conducted according to the protocols of the manufacturer by using the NEBNext R Ultra TM II DNA Library Prep Kit for Illumina R (New England Biolabs, United States). The amplification libraries were sequenced on the Illumina NovaSeq 6000 platform and 250 bp paired-end reads were generated (Guangdong Magigene Biotechnology Corporation Ltd., Guangzhou, China).

Quality Control, Sequence Assembly, and Taxonomy Assignment
Low quality ends of the pair-ended raw sequencing reads were trimmed with the assistance of fastp (Chen et al., 2018) in the "sliding and cut" way (-W = 4, -M = 20). The adapters were cut by using cutadapt (Martin, 2011). The pair-ended reads were merged by using the fastq_merge pairs function implemented in Ultra-fast sequence analysis (USEARCH) version 10 (Edgar, 2010) with the default parameters. The output merge reads, or namely "raw tags, " were undergone quality filter again by using fastp (-W = 4, -M = 20, average_qual = 25, and length_required = 200). The UNOISE3 function implemented in USEARCH was employed to denoise and filter chimera. After these steps, the optimized sequences, or "clean tags, " were clustered into operational taxonomic units (OTUs) at 97% sequence identity by the UPARSE (Edgar, 2013).
To determine valid taxonomy hits from the raw blastn output: (1) hits with identity less than 50% and coverage less than 70% of the total size of query sequence were filtered; (2) after sorting, the top 10 hits with the highest bit scores or lowest E-values were chosen to be the candidate best hits. In order to compare blast results to other tools, we introduced a statistic named "adjusted identity, " i.e., using query coverage to weight the identity. This statistics could reduce bias introduced by using only identity as cutoff (usually caused by local alignment) and provide a more intuitive statistics than E-value and bit score while comparing blast results to other tools. The adjusted identity threshold of convincing assignment by blastn was set to 97% to call a convincing assignment at genus level.
Finally, the best taxonomy assignment of all the OTUs was manually conducted based on results generated from the three tools and the inconsistence of result from these three software was manually inspected and adjusted (Supplementary Table 1).

Statistical Analysis
Rarefaction curves were generated with the online tool usearch-alpha_div_rare (version 10). 2 Microbial community compositions at the phylum level are illustrated as stacked bar charts. Hierarchical clustering of the 52 most abundant genera (OTU abundance > 0.05%) was visualized as a heatmap by using "ggplot2" version 2.2.0 (Wickham, 2016). Diversity (Shannon and Simpson indices) and richness (Chao1 and abundance-based coverage estimator (ACE) indices) estimators were calculated with 97% sequence similarity as cutoff values by the corresponding algorithms embedded in Quantitative Insights Into Microbial Ecology (QIIME) (Caporaso et al., 2010).
To illustrate the distribution pattern of microbial communities, principal component analysis (PCA; Hotelling, 1933) and principal coordinate analysis (PCoA; Zuur et al., 2007) were carried out based on Euclidean distance metrics and Bray-Curtis distance matrices, respectively. To evaluate the presence of a significant difference in the microbial community among the various sampling sites, analysis of similarities (ANOSIM; Clarke, 1993) with weighted UniFrac distance metrics or analysis of Adonis (Anderson, 2001) based on Bray-Curtis distance metrics was conducted in the vegan package embedded in R. Sparse partial least squares-discriminant analysis (sPLS-DA) was applied to reveal which organisms are responsible for the dissimilarity observed in the community composition among different samples by using the "PLS-DA" function in the R package "mixOmics" version 6.3.2 (Rohart et al., 2017).
The redundancy analysis (RDA) in combination with the Monte Carlo permutation test (ter Braak, 1990) was applied to identify a subset of environmental variables that statistically explain the greatest proportions of taxon variation. The Monte Carlo permutation test showed a p-value associated with the effect of the environmental variable on the microbial composition of the samples. This analysis was performed in R with the vegan package.
Linear discriminant analysis effect size analysis (Segata et al., 2011) was performed by using the LEfSe package embedded in an online tool 3 to identify dominant taxa among the three sites.
A phylogenetic tree of claimed OTUs with the known oil-degrading species was constructed. The best nucleotide substitution model was automatically calculated and selected with the assistance of Mega 7 (Kumar et al., 2016). Sequence alignment was conducted by muscle algorithm embedded in Mega 7. The phylogenetic tree was inferred by using the maximum likelihood method based on the Kimura 2-parameter model (Kimura, 1980).
To construct co-occurrence networks of the microbiota of different groups, pairwise intergenus correlations were calculated according to genus abundance and visualized by using the R package "igraph" version 1.2.6 (Csardi and Nepusz, 2006). The numbers of co-occurring microbe pairs under different Bray-Curtis distance matrices were compared.
Statistical procedures were carried out with R. For all the analyses, p < 0.05 was considered as statistically significant and significance levels were indicated as follows: * * * p < 0.001; * * p < 0.01; * p < 0.05.

Sampling Sites and Environmental Factors
The three sampling sites had distinct landscapes. The NS site is located outside any area with active seepage and the CH 4 concentration at this site was 0.18-0.36 µM, which was set as the baseline (Table 1). At this site, no or few faunal assemblages were observed on the seafloor ( Figure 1B). In contrast, dead clams and live tubeworms were occasionally observed at the LIS site (Figure 1C). At the HIS site, abundant live mussels and authigenic carbonates could be observed on the seafloor (Figure 1D), accompanied by continuous gas bubbling. Active seepage at the LIS and NIS sites was confirmed by the high CH 4 content in the water column measured about 50 cm above the seafloor with values of 88.3-110.4 and 1,987.1-8,831.6 µM, respectively. In addition to CH 4 concentrations, other fundamental environmental variables, including temperature, alkalinity, pH, dissolved organic carbon (DOC), nitrate, nitrite, ammonium, and phosphate concentrations, were measured; the results are listed in Table 1. Temperature, alkalinity, and pH were relatively constant among the three sites. However, DOC concentrations increased from ∼0.79 µM at the NS site to ∼1.10 µM at the LIS site and were maintained at approximately 1.50 µM at the HIS site. Interestingly, the concentrations of all the inorganic nutrients, including nitrate, nitrite, ammonium, and phosphate, were maximum at the LIS site ( Table 1). The DOC and inorganic nutrients possibly sourced from deep fluids associated with CH 4 seepage as suggested by D' Souza et al. (2016). Furthermore, the fluid seepage intensity is controlled by flow rate, thus more nutrients could dissolve in bottom water in relatively low flow rate condition at the LIS site, while nutrients upward migration to shallow water due to high flow rate at the HIS site.

Bacterial Community Composition and Alpha Diversity at the Study Sites
Deoxyribonucleic acid of microbes in the water column was extracted and sequenced. The OTU table of blank control was checked and only a few taxa could be observed with low abundance (Supplementary Table 2). Considering the extremely low counts of contaminated OTUs in the table, we concluded bacteria and archaea OTU tables with no decontamination step. The number of sequence for each sample of bacterial and archaeal libraries was included in Supplementary Table 3. An average of 53,451 and 72,919 sequences were obtained for bacteria and archaea, respectively, in all the 16 samples. Rarefaction curves based on observed bacterial and archaeal OTUs were plotted. Both the curves reached an asymptote for bacterial and archaeal libraries; therefore, the sequencing data were sufficient to detect the microorganisms present and provided reliable results (Supplementary Figure 1).
With respect to the bacterial community composition, Proteobacteria was the most abundant phylum among all the sites followed by Bacteroidetes and Actinobacteria (Figure 2 and Supplementary Table 4). The relative abundance of these phyla was significantly different between the NS and the two CH 4 seepage sites, while the LIS and HIS sites possessed similar microbial compositions with no significant difference (Supplementary Figure 2). As shown in the genus heatmap (Figure 3 and Supplementary Table 4), the distribution of top 52 abundant genera was also different between the NS and the two seepage sites, but relatively similar between the LIS and HIS sites. We also explored whether the alpha diversity of the bacterial community differed at seep sites with different seepage activities. Chao1 and ACE indices reveal community richness and estimate the OTU numbers in samples (O'Hara, 2005). Statistical analysis showed no significant differences among all the sites for the Chao1 index (p = 0.24) and ACE index (p = 0.19; Supplementary Figure 3). The Shannon and Simpson indices relate to community diversity (Schloss et al., 2011). Pairwise comparison showed that the Shannon index of the HIS site was different between the NS site (p = 0.0088), while the Simpson index was different between both the NS and LIS site (p = 0.0002; Supplementary Figure 3). Community diversity was lower when the Shannon index was smaller and the Simpson index was larger. Thus, the diversity of the HIS site was the lowest among the three sites.

Differences in Bacterial Community Structure at the Three Sites
To further explore the difference in bacterial community structure at the three sites, PCA based on Euclidean distance matrices for the genus-level compositional profiles of the three sites was performed. ANOSIM was used to assess the significance of the PCA grouping. All the sites were separated (R = 0.689, p = 0.001) along the PC1 direction (Supplementary Figure 4).
A sequential order of the NS, LIS, and HIS associated with increased seepage activity on the PC1 axis can be observed, suggesting that a dynamic pattern of bacterial community compositions might be associated with CH 4 seepage intensity. Separation was also observed in the PCoA clusters among the three sites (R 2 = 0.57, p = 0.001; Supplementary Figure 5).
Next, we identified the most discriminative taxa responsible for the difference in the bacterial community at the three sampling sites by sPLS-DA. Similar to the PCA results, the grouping of all the three sites was separated from each other (Figure 4). Genera contributing to the community variance of these sites were extracted from two components. From component 1, Alteromonas best characterized the bacterial genus compositions at the LIS and HIS sites, separating them from the NS site. Alteromonas is a prevalent marine organism able to metabolize polycyclic aromatic hydrocarbons (Jin et al., 2011). For the second component, Marinobacterium and seep sulfate-reducing bacteria-2 (SEEP-SRB2) were the characteristic genera at the HIS site. In contrast, Methylophaga, Thalassolituus, Croceicoccus, Alcanivorax, Leeuwenhoekiella, Gramella, Haemophilus, and Sphingobacterium distinguished the LIS form HIS (Figure 4). Marinobacterium is a common aerobic hydrocarbon-degrading microorganism found in oil environment (Baek et al., 2018). SEEP-SRB2 is a sulfate-reducing partner in ANME consortia for anaerobic oxidation of CH 4 (Kleindienst et al., 2012). Many Thalassolituus, Croceicoccus, and Alcanivorax species are oil degraders (Yakimov et al., 1998(Yakimov et al., , 2004Huang et al., 2015).

Environmental Variables Shaping Bacterial Community Structure
To explore the key environmental variables shaping the bacterial community composition, RDA in combination with the Monte Carlo permutation test was performed. Quantification of major geochemical parameters, including temperature, alkalinity, pH, CH 4 , DOC, nitrate, nitrite, ammonium, and phosphate, was employed and used for RDA. The Monte Carlo permutation tests indicated that CH 4 , DOC, nitrate, and ammonium had statistically significant (p < 0.05) contributions to explaining the variance in bacterial communities. These four variables together accounted for 76.99% of the explained variance. Importantly, CH 4 was the most highly significant explanatory variable (p = 0.001) and the CH 4 concentrations correlated strongly with the HIS site ( Figure 5). Therefore, CH 4 is the key environmental factor shaping the bacterial community structures at seep sites with various seepage activities.

Enriched Taxa of Study Sites With Different Seepage Activities
We applied LEfSe analysis to further determine whether specific individual bacterial taxa were differentially enriched at these sites and the results demonstrated that significantly different dominant taxa could be observed at cold seeps. Significant enrichment of genus Alcanivorax, Idiomarina, Thalassospira, Oleibacter, Marinobacter, Thalassolituus, Ruegeria, Leeuwenhoekiella, Croceicoccus, and Mesoflavibacter was FIGURE 3 | Heatmap based on the 52 most abundant genera at the three sites. The log 2 value of the taxa read number is indicated by the color gradient from light blue (low) to red (high).
identified at the LIS site, while a totally different combination of taxa, including Alteromonas, Pseudoalteromonas, Vibrio, Marinobacterium, Sulfurovum, Salinimonas, Aestuariibacter, and Thalassobius, was dominant at the HIS site (Supplementary Figure 6). Many species of the most of above genera are reported to be oil degraders, leading to the assumption of the presence of oil degraders at the two seep sites. We further confirmed this by phylogenetic analysis. We found close phylogenetic relationships appeared among known oil-degrading species with the OTUs enriched at two seepage sites belong to the genera Alcanivorax, Idiomarina, Thalassospira, Oleibacter, Marinobacter, Thalassolituus, Croceicoccus, Mesoflavibacter, Alteromonas, Pseudoalteromonas, Vibrio, Marinobacterium, Aestuariibacter, and Thalassobius. In contrast, only OTUs of four genera (Halomonas, Brevundimonas, Roseivirga, and Pseudomonas) of the 13 differentially abundant genera found at the NS site were phylogenetically closely related to oil degraders (Supplementary Figure 7). These results indicate that many putative oil degraders were enriched at the two seepage sites and natural oil seepage may occur there.

Microbial Co-occurrence Network at the Study Sites
Because bacterial communities vary from site to site, it is essential to investigate coordinated interactions of the microorganisms colonizing each CH 4 seepage site. We constructed co-occurrence networks of the microorganisms by calculating pairwise intergenus correlations based on the genus-abundance profiles for every site. Microbial co-occurrence patterns were distinct at the study sites and we calculated connections (edges) between nodes under three different Bray-Curtis distance matrices. The smaller value of Bray-Curtis distance matrices is, the stronger the co-occurrence correlation is. We found that the HIS site networks contained fewest connections between nodes under all the Bray-Curtis distance matrices. Therefore, the strength of microbial co-occurrence was weakest at the HIS site (Figure 6). FIGURE 4 | Sparse partial least squares-discriminant analysis (sPLS-DA) of the three study sites. From component 1, genera labeled in red discriminate the LIS and HIS sites from the NS site and genera labeled in blue are characteristic genera at the NS site. From component 2, genera contributed to the dissimilarity between the LIS and HIS sites are labeled with blue and red, respectively.
FIGURE 5 | Redundancy analysis (RDA) correlating several environmental variables to bacterial communities. Circles represent the bacterial community of a given sample and red arrows indicate individual environmental parameters. The concentrations of methane (CH 4 ) used for this analysis were 0.27, 99.4, and 2,207.9 µM at the NS, LIS, and HIS sites, respectively. Significance levels (marked by asterisks) were evaluated based on the Monte Carlo permutation test. ***p < 0.001; **p < 0.01; *p < 0.05.

Archaeal Community Structure at the Three Study Sites
We also analyzed archaeal taxonomic compositions at the three sites. On average, we detected 147 OTUs in the archaeal community at the three study sites; the bacterial community contained 647 OTUs, indicating that the archaeal community structure is much simpler than the bacterial community structure. Thaumarchaeota was the most abundant phylum (average 85.7%) among all the sites followed by Euryarchaeota and Nanoarchaeaeota. The relative abundance of these phyla at the three sites was similar, except for Thaumarchaeota, the abundance of which was slightly lower at the HIS site than at the NS site (Supplementary Figure 8 and Supplementary Table 4). Circos analysis was used to visualize the corresponding abundance relationships between samples and archaeal communities, showing similar relative proportions of the 15 most abundant OTUs among the three sites (Supplementary Figure 9A and Supplementary Table 5). As shown in the genus heatmap, the distribution patterns of all the genera identified were also similar among the three sites (Supplementary Figure 9B and Supplementary Table 4).
We performed further analysis of the beta diversity of the archaeal community. Clusters of three sites presented overlaps in the PCoA plot (R 2 = 0.18, p = 0.11; Supplementary Figure 10A). The PCA plot yielded similar clustering patterns with no significant separation among the three sites except for the NS and HIS sites (Supplementary Figure 10B). RDA showed that the archaeal community structure was not significantly associated with any environmental factors examined (p > 0.062; Supplementary Figure 11). Overall, the results confirmed that archaea exhibit a homogeneous distribution at different sites, indicating that they participate in general geochemical processes rather than in CH 4 metabolism in the bottom water of cold seep.

DISCUSSION
We found apart from CH 4 , which acts as the key environmental variable shaping the bacterial community structure, DOC, nitrate, and ammonium concentrations also correlated statistically significantly with bacterial community abundances. Because CH 4 oxidation by autotrophic methanotrophs leads to DOC production , the effect of DOC on shaping the bacterial community may be caused by CH 4 variation. Moreover, as many oil-degrading bacteria are enriched at the LIS and HIS sites, there is also the possibility of oil seepage at these sites with DOC produced through oil biodegradation by these bacteria (Orcutt et al., 2010). We observed highest concentrations of all the inorganic nutrients at the LIS site, which possibly sourced from deep fluids associated with CH 4 seepage (D'Souza et al., 2016). Therefore, inorganic nutrients shape the bacterial community probably due to seepage activity at the seep sites.
At the cold seeps, CH 4 oxidation likely involves communities composed of different functional guilds rather than methanotrophs alone (Dubinsky et al., 2013;Rivers et al., 2013). CH 4 levels at the cold CH 4 seeps vary intermittently from nanomolar to millimolar concentrations (Reeburgh, 2007;Ruff et al., 2013;Grupe et al., 2015). The concentration of oxygen in the bottom water will change concomitantly as CH 4 fluctuates  and oxygen availability can determine the composition of the microbial communities involved in CH 4 oxidation (Hernandez et al., 2015). We observed distinct microbial co-occurrence patterns at the study sites with the strength of microbial correlations being weakest at the HIS site. When the CH 4 concentration changes at the HIS site, the species will change and bacteria no longer present in correlation with each other; thus, the synergistic interactions between these organisms will be disrupted and results in the alteration of the topology of the co-occurrence network during CH 4 fluctuations (Shi et al., 2015). In contrast, at background site and the LIS where the CH 4 concentration did not fluctuate strongly, the strength of microbial co-occurrence was significantly greater.
Because methylotrophic bacteria may be enriched for CH 4 consumption in CH 4 seep sites, we compared canonical methanotrophs at the three sites. We found Methyloprofundus at the HIS site but not at the NS and LIS sites, while the numbers of marine methylotrophic group 2 (MMG2) and Methylotenera were similar at the three sites (Supplementary Figure 12). Members of the genus Methyloprofundus are most typically detected in sediments and within bacteriocytes of seep-associated mussels in CH 4 seeps, occasionally being found in the deepwater column (Tavormina et al., 2015). This finding indicates that Methyloprofundus may be the methanotroph responsible for CH 4 oxidation at the HIS site. However, the proportion of Methyloprofundus was relatively low at the HIS site (<0.1% of total bacterial reads). As reported by Zhang et al. (2020), methanotrophs were observed in the water of the Haima seeps but in low proportions. There is a possibility that CH 4 may be oxidized by unrecognized methanotrophs or that it may escape the CH 4 biofilter and flux to the atmosphere directly without oxidation by bacteria. Typical anaerobic cold-seep microbial groups, including methanotrophic archaea ANME-1 and ANME-2 and sulfate-reducing Desulfobulbaceae, Desulfobacteraceae, Desulfuromonadaceae, and Desulfarculaceae, were detected in extremely low proportions in our research. Similar observations were obtained in the study of Zhang et al. (2020). However, ANME-2 was found to live above a depth of approximately 600 m in the water column of the Black Sea and ANME-1 inhabited depths below 600 m; samples obtained below 600 m water depth were considered to be located in the fully anoxic water column (oxygen concentration < 1.5 µM; Schubert et al., 2006). In this study, the oxygen concentrations of bottom waters were approximately 103-109 µM (Feng et al.,in preparation 4 ). Therefore, the bottom water in the Haima seeps maintains a relatively oxic habitat that is suitable for aerobic methylotrophic bacteria for CH 4 oxidization.
Natural oil seepages have been reported in cold seeps of the Gulf of Mexico. These seepages influence the structure of the bacterial community and the abundances of many oil degraders are enhanced in oiled sediments (Orcutt et al., 2010). Several studies have demonstrated that marine obligate hydrocarbonoclastic bacteria (MOHCB) are the main actors in the degradation of hydrocarbons after oil exposure (Sanni et al., 2015). We observed many bacteria enriched at the LIS and HIS sites were putative oil degraders. Phylogenetic analysis demonstrated close phylogenetic relationship among known oil-degrading species with the OTUs enriched at the two seepage sites except for those belong to the genera Ruegeria, Leeuwenhoekiella, Sulfurovum, and Salinimonas (Yoon et al., 2012;Arahal et al., 2018;Sun et al., 2020;Tahon et al., 2020). Mesoflavibacter belongs to Bacteroidia (Okai et al., 2015), whereas Thalassospira, Croceicoccus, and Thalassobius are alphaproteobacteria (Arahal et al., 2005;Kodama et al., 2008;Huang et al., 2015); all these four genera degrade aromatic hydrocarbons. Idiomarina, Marinobacter, Alcanivorax, Oleibacter, Thalassolituus, Aestuariibacter, Alteromonas, Marinobacterium, Pseudoalteromonas, and Vibrio are gammaproteobacteria genera that consume crude oil or alkanes (Yakimov et al., 1998(Yakimov et al., , 2004Hedlund and Staley, 2001;Iwabuchi et al., 2002;Pham et al., 2009;Wang et al., 2010Wang et al., , 2014Kostka et al., 2011;Teramoto et al., 2011;Raddadi et al., 2017). Enrichment of these class was also observed by a study of the deep water horizon oil spill that occurred in the Gulf of Mexico in April 2010 with the major responders to oil intrusion in marine waters being members of gammaproteobacteria, alphaproteobacteria, and bacteroidia (Redmond and Valentine, 2012;Liu and Liu, 2013). In contrast, OTUs of only four genera (Halomonas, Brevundimonas, Roseivirga, and Pseudomonas) are capable of oil degradation at the NS site (Le Petit et al., 1975;Chaineau et al., 1999;Wang et al., 2007;Selvaratnam et al., 2016). Therefore, the observation of these putative oil-degrading bacteria in this study suggested that oil seepage may occur at cold seeps in the SCS.
We noted that oil-degrading bacteria were dominant at both the LIS and HIS sites with different seepage activities. Zhang et al. and Guan et al. studied two seep sites with different seepage intensities and showed that heterotrophic prokaryotes and UCMs were found at both the sites. These results suggest that oil may seep regardless of CH 4 seepage activity. Members of different genera occupy distinct trophic niches with diverse oil components because different MOHCB utilize different hydrocarbons; for example, Alcanivorax degrades aliphatic hydrocarbons, whereas Cycloclasticus degrades (poly-)aromatic hydrocarbons (Head et al., 2006). We noticed that the LIS and HIS sites harbored divergent hydrocarbon degraders, probably due to differences in the oil sources and components in these two seep habitats.
In addition to microbes that participate in the carbon cycle, we noted enrichment of bacteria involved in the sulfur and nitrogen cycles at both the seep sites. It has been reported that oil exposure in deep-sea sediment can also alter other metabolic pathways such as nitrogen cycling (Mason et al., 2014). It should be noted that the abundance of the denitrifier Ruegeria increased markedly at the LIS site (Supplementary Figure 6). Additionally, the thiosulfate/sulfite oxidizer Sulfurovum was dominant at the HIS site (Supplementary Figure 6), probably due to seepage of hydrogen sulfide from the sediment resulted from the sulfatedriven anaerobic oxidation of CH 4 (Lu et al., 2021).

CONCLUSION
In this study, bottom-water microbes were sampled from the NS site, LIS site, and HIS site in Haima, SCS. Our results showed that the three study sites harbor significantly divergent bacterial compositions and that the structures correlate significantly with the CH 4 concentration. The strength of the microbial correlations was weakest at the HIS site. In contrast, the archaeal community compositions were similar among the three sites and exhibited no correlation with CH 4 concentration. Diverse putative oildegrading bacteria were enriched at the two seep sites, adding new evidence that oil seepage may occur at cold seeps in the SCS. This study represents a comprehensive comparison of microbial profiles in the water column of cold seeps in the SCS, revealing that the seepage intensity has a strong impact on bacterial community dynamics.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI [accession: PRJNA715035].

AUTHOR CONTRIBUTIONS
YL and NL conceived and designed the research and wrote the manuscript. XL and ZD performed the experiments and analyzed the data. PD, JF, and JT provided the samples and associated the metadata. DC provided critical technical and scientific discussion. All authors agree with the final version of the manuscript.