Pushcore sediment samples were collected from two cold seeps (Haima: 16°43’N, 110°28′E; Site F: 22°6’N, 119°17′E) and the Xisha Trough (18°18’N, 114°08′E) located on the continent shelf of the SCS (Figure 1) during cruise TS07 by R/V “Tan Suo Yi Hao” in June 2018 as described previously (Jing et al., 2020). In total, nine stations were sampled, as follows: SQ_54, SQ_79 and SQ_81 (from Haima), SQ_62, SQ_63, and SQ_64 (from Site F), and SQ_82, SQ_84 and SQ_87 (from the Xisha Trough). The surface sediments around 0–4 cm were sliced and then immediately stored at −80°C until further analysis. In situ hydrographical parameters (i.e., temperature, depth and location) were recorded during sampling using the manned submersible, SHENHAI YONGSHI.

Chemical parameters of sediments, including total carbon (TC), total nitrogen (TN), total phosphate (TP), total sulfur (TS), nitrate, and ammonia of each site, were measured at the Institute of Mountain Hazards and Environment, Chinese Academy of Sciences (Chengdu, Sichuan, China), according to Wang et al. (2016). In total, approximately 5 g sediment was used for chemical analysis. Briefly, nitrate and ammonia were detected with a colorimetric auto-analyzer (SEAL Analytical AutoAnalyzer 3, Germany) after 2 M KCl treatment and double qualitative filter paper. TC and TN were determined by over drying the sediments at 105°C and then using an element analyzer (Elementar vario Macro cube, Germany). TS were measured by an inductively coupled plasma-optical emission spectroscopy (ICP-OES; PerkinElmer Optima 8,300, United States). TP was measured with nitric-perchloric acid using the molybdate colorimetric method with a UV2450 (Shimadzu, Japan) after digestion of the sediment (Murphy and Riley, 1962).

Triplicate samples at each station were used for DNA extraction and combined for subsequent PCR reaction and sequencing. Genomic DNA in the sediment (~0.5 g) was extracted according to the instruction manual of PowerSoil DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, United States). The extracted DNA was quantified with a NanoDrop 2000 Spectrophotometer (Thermo Scientific, Thermo Fisher Scientific, Corp.) and the quality was checked via gel electrophoresis. The DNA samples were stored at −80°C until further processing.

The V3-V4 region of bacterial 16S rRNA gene (468 bp) was amplified by PCR using the primers of 338F (5’-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′- GGACTACHVGGGTWTCTAAT-3′; Liu et al., 2016), while that of archaeal 16S rRNA gene (462 bp) was amplified with primers of 340F (5’-CCCTAYGGGGYGCASCAG-3′) and 806R (5’-GGACTACVSGGGTATCTAAT-3′; Takai and Horikoshi, 2000). Nested PCR was used for amplification of the mcrA gene (363 bp) of Methanoperedens-like archaea (ANME-2d; Vaksmaa et al., 2017), the pmoA gene (386 bp) of M. oxyfera-like bacteria (NC10 bacteria; Luesken et al., 2011), and the dsrB gene (350 bp) of Desulfosarcina-like sulfate-reducing bacteria (DSRB; Rampinelli et al., 2008). The information of specific PCR primers for above functional genes are listed in Table S1. PCR products were examined with SYBR Safe stained 1.2% agarose gels. The paired-end sequencing of the all amplicons were then performed with an Illumina HiSeq PE250 sequencer (Novogene Co., Ltd.).¹

To obtain high-quality sequencing data and improve the accuracy of subsequent bioinformation analysis, adapters at both ends of the sequences were removed firstly using the q2-cutadapt plugin after demultiplex; Then, the DADA2 (v1.16) plug-in of QIIME2 v2020.2 was used to filter, dereplicate, identify chimeric sequences, and merge Paired-end (PE) reads (Callahan et al., 2016). PE reads were merged with a minimum length of 12 bp overlap and the representative sequences were picked. The unique amplicon sequence variants (ASVs) table was filtered out by q2-filterfeature after removing the ASVs with frequencies less than 10. The resulting representative ASV sequences for functional genes (mcrA, pmoA and dsrB) were used to calculate phylogenetic trees for subsequent analyses.

Taxonomic classification of 16S rRNA (bacteria and archaea) was processed using the q2-classifysklearn algorithm, and the SILVA (V.132) database was used as a reference with a threshold of 0.8. Annotations were obtained after removing contamination using the q2-feature-table plugin and visualized by the q2-taxa-barplot plugin. The ASVs annotated as mitochondria, chloroplasts, or eukaryotes were also removed using the qiime taxa filter-table and qiime taxa filter-seqs plugins of QIIME2. To detect potential biomarkers, linear discriminant analysis (LDA) effect size (LEfSe) statistical analysis was performed on the Galaxy platform (Segata et al., 2011).² For the prediction of functional and metabolic profiles of the bacterial and archaeal community based on the 16S rRNA gene sequences, the recently developed open-source R package Tax4Fun (Asshauer et al., 2015) was used with the short reads mode disabled along with the SILVA database 123 as required.

Alpha diversity (Shannon and Chao1 indices) and beta diversity were generated with QIIME 2 using q2-diversity, then visualized using box plots and non-metric multidimensional scaling (nMDS) plots at the class level. Pearson’s correlation coefficients were calculated to identify a possible differentiation of the communities under the constraint of environmental factors, and assess correlations between environmental variables and community variability. ANOSIM (analysis of similarities) was used to analyze the similarities of the microbial community compositions among different regions at the class level.

A phylogenetic analysis of functional gene (mcrA, pmoA and dsrB) sequences was performed using Mega X (Kumar et al., 2018). Reference sequences were retrieved from GenBank database using nucleotide tool BLAST, Popset, and Batch Entrez at NCBI.³ The gene sequences were aligned using the ClustalW algorithm of Mega X (Kumar et al., 2018). After testing the best substitution model, phylogenetic trees were constructed using the maximum likelihood method with the substitution model (Tamura-Nei model for mcrA, Tamura 3-parameter model for pmoA and dsrB) with a bootstrap value of 1,000. The phylogenetic tree was visualized and edited through ITOL (V2.0; Letunic and Bork, 2006).

The abundance of the functional genes (mcrA, pmoA and dsrB) was quantified using the Bio-Rad System (Bio-Rad Inc., United States) and TB Green Premix^® Ex Taq II (Takara Bio Inc., Shiga, Japan) with primers McrA159F/McrA345R (Vaksmaa et al., 2017), cmo182F/cmo568R (Luesken et al., 2011) and DSRp2060F/DSR4R (Geets et al., 2006) respectively. The specific information for all the qPCR primers were listed in Table S1. Standard curves were constructed using a series of tenfold dilutions of the standard plasmids (known copy number) containing the targeted genes. Triplicate qPCR reactions were performed for each sample with double-distilled water as a negative control, and the gene copy number was normalized to the quantity of the gene.

The Xisha trough was located in between the two seep regions, and geographically closer to Haima (Figure 1). The highest contents of TN (0.09% ~ 0.17%) and TC (0.67% ~ 1.27%) were found in Haima (Table 1). The carbon to nitrogen ratios (C/N; 6.87 ~ 8.30) were higher in the cold seeps than in the trough, and the highest values were found at Stn. SQ64 in Site F. The concentration of NH₄⁺ (9.49 ~ 17.10 mg/kg) was highest at Stn. SQ64 in the Site F and lowest at Stn. SQ81 in the Haima. The highest and lowest concentration of NO₃⁻ (1.01 ~ 1.41 mg/kg) was detected at respective Stns. SQ87 and SQ82 in the Xisha trough. TP (477.01 ~ 1,071.21 mg/kg) was significantly higher in Site F, especially at Stn. SQ64. In general, Site F contained high concentrations of TP, NH₄⁺ and C/N ratio, whereas Haima had high TN and TC contents.

The sequencing process yielded 460,109 quality reads for bacterial 16S rRNA gene and 549,475 quality reads for archaeal 16S rRNA gene from 9 stations (Table S2). A total of 9,251 ASVs from bacterial 16S rRNA gene and 5,234 ASVs from archaeal 16S rRNA gene were identified in all nine samples. The bacterial ASVs were classified into 18 phyla and 34 classes, and the archaeal ASVs were classified into 9 phyla and 9 classes (Figures 2A,B). The average alpha diversity (Shannon) of bacterial community at the class level was significantly higher than that of archaeal community in the three regions, especially in the Site F and the Xisha trough (p < 0.01, Figure 2C). NMDS analysis of bacteria and archaea at the class level showed that three distinct clusters were formed, corresponding to the three sampling regions (Figures 2D,E).

The major bacterial phyla in all samples were Proteobacteria, Chloroflexi and Acidobacteria (Figure 2A). Among them, Proteobacteria was the predominant bacterial phylum, and comprised mainly of α- and γ-Proteobacteria (Figure 2A). In addition, NC10 bacteria belonging to Rokubacteria were detected in all study areas (Figure 2A). Within the 9 archaeal phyla, Thaumarchaeota was the most abundant archaeal phylum in all samples, and dominated by Nitrososphaeria (Figure 2B). Thermoplasmata and Methanomicrobia belonging to Euryarchaeota, together with Bathyarchaeia (Crenarchaeota) occupied a higher proportion in Haima especially at Stn. SQ81 (Figure 2B). The relative abundance of Marine Benthic Group A, Thaumarchaeota and Thermoplasmata, was higher in the trough than the cold seeps, while that of Methanocellales and Methanomicrobia was higher the in cold seeps than the trough (Figure 2B).

LEfSe analysis demonstrated that NC10 bacteria (Rokubacteria) and Methylococcales accounted for higher proportion in the trough and the cold seeps, respectively (p < 0.05, Figure 3A; Table S3). The relative abundance of Marine Benthic Group A, Thaumarchaeota and Thermoplasmata, was higher in the trough than in the cold seeps, while that of Methanocellales and Methanosarcinales showed the opposite trend (p < 0.05, Figure 3B).

For functional prediction of bacteria with Tax4fun based on total ASVs of bacterial 16S rRNA gene, the abundance related to acetoclastic methanogens was the highest in all methanogenesis pathways, especially in the Xisha trough, while that related to CO₂ and methylotrophic methanogens was higher in the cold seeps (Figure S1A). The genes related to denitrification and dissimilatory nitrate reduction was more abundant in the Xisha trough, while that of assimilatory sulfate reduction was more abundant in Site F (Figure S1A).

For functional prediction of archaea with Tax4fun based on total ASVs of archaeal 16S rRNA gene, the abundance of methane metabolism was the highest in the cold seeps especially in Haima, while that of denitrification, dissimilatory nitrate reduction and assimilatory sulfate reduction was higher in the Xisha trough (Figure S1B).

A total of 638,497 high-quality mcrA gene sequences were obtained from all samples excluding Stn. SQ81, and were classified into 222 ASVs using DADA2 (Table S2). Top 53 ASVs (accounting for >95% of the total retrieved mcrA gene sequences) were used to construct a phylogenetic tree, and fell into five distinct clusters, i.e., ANME-2e, ANME-2c, ANME-2d, Methanosarcinaceae, and Methanocellales (Figure 4A). ANME-2e (22ASVs, 27.9%) and ANME-2d (18ASVs, 27.6%) were the two major clusters.

A total of 624,726 high-quality pmoA gene sequences were obtained in all nine samples, and were classified into 251 ASVs using DADA2 (Table S2). Top 20 ASVs (accounting for >85% of the total retrieved pmoA gene sequences) were used to construct a phylogenetic tree, with two distinct clusters formed (Figure 4B). Cluster I contained most of the ASVs (15 of 20 ASVs).

A total of 696,892 high-quality dsrB gene sequences were obtained in all nine samples and were classified into 2,646 ASVs using DADA2 (Table S2). Seven distinct clusters were formed based on the top 20 dsrB ASVs with affiliated sequences from the GeneBank database, i.e., Group I, Syntrophobacteraceae, Desulfobacteraceae, Group II, Desulfobulbaceae, Group III and Group IV (Figure 4C).

For mcrA gene, the highest proportion of Methanosarcinaceae and lowest proportion of ANME-2e were detected in Site F (Figure 4A). ANME-2c and Methanocellales were major group in the Haima and the Xisha trough, respectively. The relative abundance of ANME-2d was higher in the cold seeps (accounting for 23.3% in the Haima and 42.4% in the Site F, respectively) than that in the trough (accounting for 4.1%; Figure 4A). Meanwhile, distinct ANME-2d predominant in each sample, e.g., ASV5, ASV4 and ASV24 was mainly distributed in the Haima, the Site F and the Xisha trough, respectively, (Figure S1C). For pmoA gene, Cluster I affiliated with NC10 bacteria accounted for higher proportion in the cold seeps than the trough. ASV1 dominant in all samples, while ASV2 was mainly in the cold seeps (Figure S1D). For dsrB gene, Group I and Syntrophobacteraceae were mainly in Site F (respective of 44.3 and 36.5%); higher relative abundance of Desulfobacteraceae, Group III and Group IV were found in the Xisha trough, while Group II and Desulfobulbaceae accounted for higher proportion in the Haima (Figure 4C; Figure S1E).

In terms of alpha diversity, there was no significant difference of mcrA group between the cold seeps and the trough (Figure 5A), while significantly higher diversity of pmoA and dsrB groups was present in the trough and the cold seeps, respectively (p < 0.05, Figures 5B,C). Regarding gene abundance quantified by qPCR, significantly higher abundance of mcrA, pmoA and dsrB gene was detected from the cold seeps than from the trough (p < 0.05, Figure 5D). NMDS analysis of functional genes showed that three distinct clusters corresponding to the different sampling regions were formed (Figure 5E).

Pearson’s correlation coefficients demonstrated that significant correlation existed for Group II/III of dsrB with TS and C/N Ratio, and for ANME-2c/2d of mcrA with TS and TP (p < 0.05, Table 2). Regarding alpha diversity, it was significantly correlated with TN and TC in the bacterial community, and with TN, TC and TS in the archaeal community (p < 0.05, Table 2).

Bacterial community exhibited higher diversity than that of archaeal community in the three studied regions in this study, consisting with the finding of previous studies conducted in the cold seeps of the SCS (Zhang et al., 2012; Cui et al., 2019; Jing et al., 2020). Compared with trough, higher proportions of SRB (Desulfobulbus, Desulfococcus and Desulfobacteraceae) were found in the cold seeps, where typically syntrophic consortium formed by those SRB groups with anaerobic methane oxidizers (Boetius et al., 2000; Orphan et al., 2001). On the other hand, NC10 bacteria accounted for higher proportion in the trough, where has been proposed to support a widespread occurrence of this bacterial group (Chen et al., 2014) with the potential of performing N-DAMO pathway with the co-operation of anammox bacteria (Hu et al., 2009; Ettwig et al., 2010). In addition, Methylococcales (type I methanotrophs) occupied a high proportion in the cold seeps, very possibly be favored by relatively high methane released after AOM process (Hanson and Hanson, 1996).

As for archaea, the prevalence of methylotrophic Methanosarcinales and acetogenotrophic Methanocellales in the cold seep sediments suggested a potentially high level of in situ methane production, and in agreement with the functional prediction based on 16S rRNA genes in this study. In addition, archaeal groups, i.e., Methanoperedenaceae (ANME-2d), ANME-3, ANME-2a/2b and ANME-2c, were mainly detected in the cold seeps. It was already known that ANME −2a, −2b, −2c, and − 3 could form a consortium with SRB to participate in the SAMO process (Knittel and Boetius, 2009), therefore, their presence might further prove that SAMO processes widely occur in the hydrate-bearing seeps of the northern SCS (Niu et al., 2017; Cui et al., 2019). The detection of ANME-2d clades Methanoperedens in the cold seeps indicated that the DAMO-related microbes existed in the hydrate-bearing cold seeps (Haroon et al., 2013), although their real ecological role in this process still need direct proofs. Different ANME clades had varied spatial distribution, for example, ANME-3 predominant in this study, while ANME-1 and ANME-1b predominant in the respective of Jiulong methane reef area (Cui et al., 2016) and the GMGS2 drilling area of the SCS (Cui et al., 2019). This reflected the complexity and heterogeneity of the microbial communities in different cold seeps, and the necessity of conducting specific studies on each individual seep.

Based on the phylogeny of mcrA gene ANME-2d cluster affiliated with sequences recovered from Yangtze Estuary sediment (Zheng et al., 2020), wetland (Chen et al., 2017) and Indonesian river sediment (Vaksmaa et al., 2017), which were typical habitats for the occurrence of methanogenesis (Wang et al., 2015) and Nr-DAMO processes (Haroon et al., 2013). Higher proportion of ANME-2d in the cold seeps was consistent with the findings based on archaeal 16S rRNA gene analysis in this study, although this archaeal cluster was undetectable in the Haima cold seep previously using another conventional mcrA primer (Niu et al., 2017). As for pmoA gene, two distinct clusters formed, with Cluster I was closely related to sequences recovered from the mudflat in Mai Po Nature Reserve of Hong Kong (Chen et al., 2015a) and sediments of the SCS (Chen et al., 2014, 2015b). NC10 bacteria have been reported using pmoA gene from the sediments in the SCS previously (Chen et al., 2014, 2015b, 2022). The association of ANME-2d archaea and NC10 bacteria with DAMO process has been proposed (Chen et al., 2022), future studies with isotopic tracing experiments will be helpful to elucidate their functions and relative contributions to this process.

Comparatively, more distinct clusters were formed for dsrB gene, suggesting highly diversified dsrB genotypes existed. Among the seven clusters, clusters of Syntrophobacteraceae, Desulfobulbaceae and Group IV were closely related to sequences from deep-sea sediments of Nankai Trough (Kaneko et al., 2007); Desulfobacteraceae cluster was closely related to sequences from sediments of Changjiang Estuary (He et al., 2015); clusters of Group II/III were closely related to sequences from coastal marine sediment of Aarhus Bay (Petro et al., 2019). The high diversity of SRB genotypes might be associated with the different ANME groups detected in the same region. This phenomenon has also been observed in the pockmark (Lazar et al., 2011) and Thuwal cold seeps (Lee et al., 2014), might be influenced by the sulfate concentrations and organic substrates availability (Leloup et al., 2007; Lazar et al., 2011). Highly diversified substrate utilizing capability of different SRB groups reflected a better adaptability to the deep-sea environment (Dhillon et al., 2003), and a potentially higher ecological contribution to AOM progress in the cold seeps and trough.

Higher gene abundance of all the three functional groups in the cold seeps than in the trough might suggest that those functional groups were selectively favored by the enriched substrates in the cold seeps (Zhou et al., 2015; Niu et al., 2017; Cui et al., 2019). In the Xisha trough, significantly higher abundance of NC10 group might be associated with the frequently predicted functions of denitrification and dissimilatory nitrate reduction. This suggested that N-DAMO process was more important in the hydrate-bearing trough, although its contribution to the AOM process should be estimated using the isotopic tracing and active microbial assemblages involved deserve further investigation at the RNA level.

Distinct prokaryotic and functional microbial communities were formed in the three respective hydrate-bearing locations. Distinct microbial communities formed in different cold seeps have been reported worldwide, and were locally selected by the biotic and abiotic factors (Ruff et al., 2015; Vigneron et al., 2019). The degree of microbial endemism in the methane seep suggests a high local diversification in the heterogeneous cold-seep ecosystems (Ruff et al., 2015). The clear biogeographical distribution pattern of microbial communities were attributed to the in situ geochemical conditions in the sediments, especially TN, TC, TP, NO₃⁻, and TS; and their importance driving the spatial distribution of microbial communities in the cold-seep ecosystems have been reported previously (Heijs et al., 2007; Roalkvam et al., 2011; Shao et al., 2014; Niu et al., 2017). However, by far different environment parameters were measured in different studies, and it is difficult to ascertain which one was the most important factor. Previous studies indicated strong correlations of microbial communities in the cold seep sediment with the concentrations (Zhang et al., 2012) and physical forms of methane (Cui et al., 2019). In those hydrate-bearing ecosystems, it would be reasonable to assume that the concentrations and physical forms of methane as the key drivers, although no in situ methane concentration was reported in most related studies. Methane concentrations were not measured in our study, but bubbling of fluids from the cold seeps were observed during diving with the deep-sea HOV in the cruise. Haima has exhibited a decline in activity in recent years (Liang et al., 2017), while Site F is currently active (Feng and Chen, 2015) with high concentration of methane gas (Lin et al., 2007), and the Xisha trough with large amount of gas hydrate, without cold seeps formed currently (Zhang et al., 2019). Therefore, concentrations and physical forms of methane in those different regions would be different, subsequently supporting distinct prokaryotic and functional communities, and should be measured and recorded in the future studies.

16S rRNA genes as phylogenetic markers for amplicon sequencing have become routine in the microbial ecology studies in recent years (Woese and Fox, 1977), but functional genes that are indicative of a particular microbial guild were applied increasingly in recent years (Zeleke et al., 2013). Both 16S rRNA gene and functional genes were applied in this study. It seems 16S rRNA gene was more suitable for routine community composition analysis, however, to get a better resolution of functional groups, the specific primers targeting on the functional genes would be necessary. For example, different ANME clades (referred to as ANME −2a, −2b, −2c, and − 3) were revealed by the mcrA gene, while those clades together occupied a relatively lower proportion in the whole archaeal community (ranging from 0.04 to 1.70%) based on the archaeal 16S rRNA gene (Table S3). In addition, the relative abundance of SRB was extremely low in all samples based on the analysis of bacterial 16S rRNA genes. In fact, SRB was known belonging to several diverse clusters, a combination of different primers or probes are needed to cover entire SRB communities specifically (Lücker et al., 2007), and 16S rRNA gene-based monitoring of SRB is particularly not sufficient. Therefore, it is necessary to use specific primers of functional genes combined with 16S rRNA to get a comprehensive picture of the microbial communities, and this is particularly true for the environmental samples containing complex microbial gene pools (Lücker et al., 2007; He et al., 2015).