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₄ 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).¹ 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.

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 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^® 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₂O) to a total of 50 μL. The “contaminant” control was performed with sterile ddH₂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.^® Gel Extraction Kit (Omega, United States). Library construction was conducted according to the protocols of the manufacturer by using the NEBNext^® Ultra^TM II DNA Library Prep Kit for Illumina^® (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).

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).

Taxonomy assignment was conducted by using a combination of three tools: (1) blastn (word_size = 21) embedded in standalone version of Basic Local Alignment Search Tool (BLAST+) version 2.11.0 (United States) (Camacho et al., 2009) against the SILVA 16S/18S (SSU) reference database version 138.1 (Bremen, Germany) (Quast et al., 2013); (2) OutTab function implemented in the USEARCH (identity cutoff = 97%) against its embedded Ribosomal Database Project (RDP) taxonomy 18 database; and (3) and the RDP classifier version 2.13 (confidence cutoff = 0.8; Wang et al., 2007) against SILVA SSU reference database version 138.1.

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).

Rarefaction curves were generated with the online tool usearch-alpha_div_rare (version 10).² 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³ 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.

The three sampling sites had distinct landscapes. The NS site is located outside any area with active seepage and the CH₄ 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₄ 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₄ 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₄ 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.

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₄ 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.

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₄ seepage intensity. Separation was also observed in the PCoA clusters among the three sites (R² = 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₄ (Kleindienst et al., 2012). Many Thalassolituus, Croceicoccus, and Alcanivorax species are oil degraders (Yakimov et al., 1998, 2004; Huang et al., 2015).

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₄, DOC, nitrate, nitrite, ammonium, and phosphate, was employed and used for RDA. The Monte Carlo permutation tests indicated that CH₄, 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₄ was the most highly significant explanatory variable (p = 0.001) and the CH₄ concentrations correlated strongly with the HIS site (Figure 5). Therefore, CH₄ is the key environmental factor shaping the bacterial community structures at seep sites with various 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 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.

Because bacterial communities vary from site to site, it is essential to investigate coordinated interactions of the microorganisms colonizing each CH₄ 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).

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² = 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₄ metabolism in the bottom water of cold seep.