Sediment push cores were obtained from four habitats of the Haima cold seep in the north of SCS (16.9°N, 110.4°E) via remotely operated vehicle (ROV) in September 2020. The non-seepage habitat (NS) was normal seafloor (1,393 m, Figure 1A), representing the status before seepage. Clam bed (CB) inhabited Vesicomyidae clams, and microbial mats (1,362 m, Figure 1A). MB and MV were both mussel habitats. MB possessed large mussel aggregates on muddy sediment (1,394 m, Figure 1A). Gas venting and large carbonate crusts were found in the dominant mussel MV (1,440 m, Figure 1A). One push core was retrieved from NS (25 cm), CB (35 cm), and MV (24 cm), respectively, and two push cores were obtained from MB (24 and 23 cm). Porewater was extracted using a syringe under vacuum upon retrieval at 4°C, preserved with a 10 μL saturated HgCl₂ solution, and stored at 4°C under an N₂ headspace. Then, the push cores were cut into 3 cm per sample and stored at −80°C before use. The methane concentration of the bottom seawater was recorded using a sensor on the Haima ROV. The maps were illustrated with Ocean Data View version 5.4.0 (Schlitzer, 2015).

Geochemical variables were analyzed following our previous study (Luo et al., 2015). Briefly, porewater samples were diluted 1:500 and 1:200 with ddH₂O to determine the concentrations of negative ions [sulfate (SO42-) and chloride (Cl⁻)] and positive ions [magnesium ion (Mg²⁺) and calcium ion (Ca²⁺)], using a Dionex ICS-900 ion chromatography. An Ion Pac AS23-type column and a 4.5 mM Na₂CO₃ and 0.8 mM NaHCO₃ mixed solution in an anion system were used to determine SO42- and Cl⁻ concentrations. An IonPac CS12A-type column and an 11 mM H₂SO₄ solution in a cation system were used to determine magnesium [Mg²⁺] and calcium [Ca²⁺] concentrations. Concentrations of nitrate [NO3-], nitrite [NO2-], ammonia [NH4+], silicate [SiO32-], and phosphate [PO42-] were measured using a Lachat QuickChem8500 autoanalyzer (Lachat Instruments, Loveland, CO, USA) following standard colorimetric methods (Grasshoff et al., 2009).

Total DNA was extracted from 2 g of sediments using the TGuide S96 Magnetic Soil/Stool DNA Kit [Tiangen Biotech (Beijing) Co., Ltd.] for full-length 16S rRNA gene amplicon and metagenomic sequencing according to the manufacturer's instructions. The DNA concentration of the samples was measured with the Qubit dsDNA HS Assay Kit and Qubit 4.0 Fluorometer (Invitrogen, Thermo Fisher Scientific, Oregon, USA).

The bacterial and archaeal 16S rRNA genes were amplified using the 27F/1492R primer pairs (Delong, 1992) and 20F/1492R primer pairs (Delong, 1992; Chen et al., 2016), respectively. PCR progress was shown as follows: 95°C for 2 min; followed by 25 cycles of 98°C for 10 s, 55°C for 30 s, and 72°C for 90 s; and finished with 72°C for 2 min. Amplicons were purified with Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN). SMRTbell libraries were prepared with the SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences) and sequenced on the PacBio Sequel II platform. Raw data were primarily filtered and demultiplexed using the SMRT Link software (version 8.0) with minPasses ≥5 and minPredictedAccuracy ≥0.9 to obtain the circular consensus sequencing (CCS) reads. The LIMA (version 1.7.0) was employed to assign the CCS to the corresponding samples based on the barcodes. CCS reads without primers and with a length over the range of 1,200–1,650 bp were discarded, and quality filtering was performed using Cutadapt (Martin, 2011; version 2.7). Amplicon sequence variants (ASV) were generated with DADA2 (Callahan et al., 2016) in QIIME2 (Bolyen et al., 2019), and the ones with an abundance of < 0.005% were filtered (Bokulich et al., 2013). The ASV matrices were rarefied by the minimum sample size for archaea (17,388) and bacteria (6351). Taxonomy annotation was performed based on the Naive Bayes classifier in QIIME2 (Bolyen et al., 2019) using the SILVA database (Quast et al., 2012; release 132).

Alpha diversity indices were calculated with the picante (Kembel et al., 2010) package in R. Bray-Curtis Similarity Index-based non-metric multidimensional scaling (NMDS) was calculated via the vegan package (Dixon, 2003) in R. Analysis of similarity (ANOSIM), permutational multivariate analysis of variance (PERMANOVA), Mantel test, and Student's t-test were performed via vegan (Dixon, 2003) package (with 9999 permutations for ANOSIM, PERMANOVA, and Mantel test). Distance-based redundancy analysis was performed via a vegan package. Spearman correlation was calculated with the psych (Revelle, 2011) package in R. For phylogenic analysis, the sequence of 16S rRNA gene ASVs was aligned using the MUSCLE algorithm (Edgar, 2004). Maximum likelihood trees were constructed using the Tamura-Nei model with all sites, bootstrapped with 500 replicates by MegaX (Kumar et al., 2018).

To identify the interactions of microbial species at different aggregates, molecular ecological networks (MENs) were constructed on the basis of the Pearson correlations of bacterial and archaeal ASV abundances, followed by an RMT-based approach that determines the correlation cut-off threshold in an automatic fashion (Luo et al., 2006, 2007; Zhou et al., 2010) with the Molecular Ecological Network Analyses Pipeline (Deng et al., 2012). Archaeal and bacterial ASVs from the same sample were pooled for the construction of MENs (6–9 samples for each man, Supplementary Table 1), and only the ASVs present in at least 6 samples were included for correlation calculation (Deng et al., 2012).

The topological indices, including total nodes, total links, power-law fitting of node degrees, average K (avgK), average CC (avgCC), average path distance (geodesic distance, GD), geodesic efficiency (E), modularity, harmonic geodesic distance (HD), centralization of stress centrality (CS), centralization of betweenness (CB), centralization of eigenvector centrality (CE), density, modularity, and efficiency, were calculated in the MENAP interface (Deng et al., 2012). Random networks were generated for each empirical network by randomly rewiring the links among the nodes while constraining nodes and links, following the Maslov–Sneppen procedure (Maslov and Sneppen, 2002) in the MENAP.

A total of 10 nanograms of DNA extracted from 2 g of sediments using the TGuide S96 Magnetic Soil/Stool DNA Kit [Tiangen Biotech (Beijing) Co., Ltd.] was used to produce a library via the VAHTS^® Universal Plus DNA Library Pren Kit for Illumina, and the pooled libraries were sequenced on the Illumina Novaseq 6000 platform. Raw data were primarily filtered using Trimmomatic (Bolger et al., 2014; version 0.33) with the following parameters: Pe leading 3, trailing 3, sliding window 50:20, and minLen 120. Filtered reads were assembled via MEGAHIT (Li et al., 2015; version v1.1.2) for each sample, and short contigs (< 300 bp) were filtered. Identification of the opening reading frame (ORF) was conducted using MetaGeneMark (http://exon.gatech.edu/meta_gmhmmp.cgi, version 3.26) with the default parameter. Non-redundant genes were clustered in line with a sequence similarity of 95% and a coverage of 90% using MMseq2 (https://github.com/soedinglab/mmseqs2, version 11-e1a1c). Non-redundant genes were annotated with the Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-Prot, and NR databases in DIAMOND (version 0.9.29, cutoff: e-value 1e-5). To construct the metabolic pathway, genes with an average relative abundance of over 1e-06 were mapped via the Reconstruct tool in KEGG (Kanehisa et al., 2021). The schematic diagram was generated via BioRender (https://app.biorender.com) and modified manually.

A total of four aggregates, namely, non-seepage muddy seafloor (NS), clam bed (CB), mussel bed (MB), and mussel aggregates with methane vent and carbonate rocks (MV), were found in Haima's cold seep during a cruise in 2020. Dead clams were found in the mussel habitats (MB and MV), but mussels (either living individuals or remnants) were not found in CB (Figure 1A). Accordingly, based on our observations and previous studies (Bergquist et al., 2003; Cordes et al., 2006; Bowden et al., 2013), the aggregates may represent different stages of cold seep development. Geochemical variables in the sediment showed evident differences among aggregates (Figure 1B, Supplementary Figure 1). Sulfate and calcium concentrations showed dramatic vertical declines in the push cores obtained from MB and MV (Supplementary Figure 1). The methane concentrations for bottom seawater were recorded as 0, 1.1, 5.6, and over 9.9 μmol for NS, CB, MB, and MV, respectively, indicating an increase in methane flux from NS to MV. Ammonium and phosphate, usually derived by methanogenesis (Quevedo et al., 1996; Lippens and De Vrieze, 2019), increased in concentration from NS to MV. The decreased concentrations of sulfate and calcium from NS to MV indicated an increase in AOM intensity. The concentration of silicate increased from CB to MV, indicating a more alkaline environment that promoted its solvation (Van Cappellen and Qiu, 1997; Smrzka et al., 2015). Variation of geochemical variables was not found in a vertical direction, except in MB. The dramatic decline of sulfate in the subsurface samples indicated that they were in the sulfate methane transition zone (SMTZ). The geochemical metadata exhibited a significant difference between aggregates and an increase in heterogeneity in MB.B and MV (ANOSIM, R = 0.61, p < 0.001, Figure 1C). Therefore, we divided the sediment samples into five groups (NS, CB, MB.R, MB.B, and MV) according to the geochemical variables (Supplementary Table 1).

A total of 1,099,280 full-length bacterial 16S rRNA gene sequences from 40 samples and 676,132 archaeal sequences from 38 samples were retrieved (Supplementary Table 1). The Shannon index of the archaeal community decreased from NS to MV (Figure 2A), while that for the bacterial community increased from NS to MB but decreased at MV (Figure 2B). Non-metric multidimensional scaling (NMDS) based on Bray–Curtis distance revealed a clear alternation in archaeal and bacterial community structure in different aggregates of Haima's cold seep (Figures 2C, D), which was also supported by PERMANOVA and ANOSIM (Supplementary Tables 3, 4). Anaerobic methanotrophic archaea (ANME) affiliated with Euryarchaeota were the predominant archaea in CB, MB, and MV, comprising 81.3–98.4% of the archaeal community. The ANME-2b clade was hardly found in NS, MB, and MV but dominant in MV (Figure 2E). For bacterial communities, the abundance of Proteobacteria, Atribacteria, and Campylobacterota increased from 28.7% in NS up to 72.7–89.6% in CB, MB, and MV. Proteobacteria were the most diverse bacterial phylum and mainly consisted of sulfate-reducing bacteria (SRB). Syntrophobacteraceae were dominant and prevalent in NS, while SEEP-SRB1, SEEP-SRB2, SEEP-SRB4, Desulfocapsa, and Desulfatiglans were dominant in CB and MB. Milano-WF1B-44 became the predominant Proteobacteria in MV, whose relative abundance rose to 38.4%. JS-1 clade affiliated with Atribacteria largely inhabited CB and MB, comprising 9.1–65.5% of bacteria. The sulfide-oxidating bacteria (SOB) Sulfurovum (Assié et al., 2020), affiliated with Campylobacterota, has been abundant since CB, comprising 5.3–61.0% of bacteria (Figure 2F).

The community diversity and its correlation with environmental variables were evaluated for the four aggregates. Alpha diversity indices of the archaeal community were largely correlated with concentrations of sulfate, calcium, ammonium, and silicate (Spearman, |R| > 0.3, p < 0.05, Supplementary Table 5). However, the diversity indices of the bacterial community were correlated with concentrations of nitrite and nitrate (Spearman, |R| > 0.3, p < 0.05, Supplementary Table 6). The composition of the archaeal community was strongly correlated with concentrations of sulfate and silicate (Mantel test, R > 0.2, p < 0.05) and exhibited a weak but significant correlation with magnesium (Mantel test, 0.1 < R < 0.2, p < 0.05). In contrast, the composition of the bacterial community was strongly correlated with concentrations of sulfate and calcium (Mantel test, R > 0.2, p < 0.05) and exhibited a weak but significant correlation with concentrations of magnesium, ammonium, nitrate, and nitrite (Mantel test, 0.1 < R < 0.2, p < 0.05, Figure 2G).

To evaluate the associations among microbial species, molecular ecological networks (MENs) were constructed using a random matrix theory (RMT)-based approach (Luo et al., 2006, 2007; Zhou et al., 2010). The network nodes represent microbial species, and links between nodes represent their abundances that were significantly correlated. All the networks displayed scale-free, small-world, and modular properties. They exhibited significant differences in topological indices from the corresponding random networks (Supplementary Table 7), suggesting that the constructed networks were ecologically meaningful. We found the number of network nodes and modules that peaked in MB.R (Figure 3A, Supplementary Table 7), which was in line with changes in bacterial species diversity (Figure 2D). However, the network nodes were most densely correlated in MV, indicating that the species were highly associated in this aggregate. The proportion of negative links was the lowest at CB but gradually increased to 52.8% in MV (Figure 3B). Such increases in negative link proportion were also found among ANME and SRB species, as well as between ANME, SRB, JS-1, and Sulfurovum species (Supplementary Figure 2). Moreover, topological indices of the networks, including the proportion of negative links (Neg), average degree (avgK), geodesic efficiency (E), and density, were positively correlated with silicate, while modularity, average path distance (GD), harmonic geodesic distance (HD), centralization of betweenness (CB), centralization of eigenvector centrality (CE), and efficiency were negatively correlated with silicate (Spearman, |R| > 0.8, p < 0.05, Figure 3C).

To evaluate the distribution patterns and their correlations with geochemical variables, phylogenetic trees were constructed for networked species affiliated with ANME, SRB, JS-1, Sulfurovum, and Milano-WF1B-44. ANME-1b was the dominant networked archaea, and most of them inhabited CB. Interestingly, the abundance of ARC1 kept increasing and peaked at MB. ARC31 and ARC111, affiliated with ANME-2b, were prevalent in MV and one of the dominant species. Other ANME species did not show an apparent aggregate-specific distribution pattern. The distribution of ANME species was largely positively correlated with sulfate (Spearman, R > 0.3, p < 0.05) but negatively correlated with silicate (Spearman, R < −0.3, p < 0.05). In contrast, ARC1, ARC31, and ARC111 have shown a negative correlation with sulfate and calcium (Spearman, R < −0.3, p < 0.05) but a positive correlation with silicate and phosphate (Spearman, R > 0.3, p < 0.05), indicating their prevalence in methane and correlation with AOM (Figure 4A). The networked JS-1 species, ubiquitous cold seep taxa involved in hydrocarbon fermentation (Liu et al., 2019; Dong et al., 2020), were specific in the CB and MB. Their distribution was dominantly correlated with sulfate and silicate (Spearman, |R| > 0.3, p < 0.05), indicating their sensitivity to the geochemical changes caused by AOM (Supplementary Figure 3).

SRB directly participates in the biogeochemical cycle of carbon and sulfur (Fukui et al., 1999; Kleindienst et al., 2014; Jaekel et al., 2015; Davidova et al., 2019), showing a clear aggregate-specific distribution pattern in Haima's cold seep. Syntrophobacteraceae, mediating the transformation of long-chain alkane to carbon dioxide or acetate (Tan et al., 2014a,b), were the dominant SRB in NS but were absent in MB.B and MV. SEEP-SRB1 and SEEP-SRB2, the well-known ANME partners (Kleindienst et al., 2012), were prevalent in CB and MB. Other networked SRB species, such as Desulfocapsa, Desulfatiglans, Desulfoconvexum, and Desulfobulbus, were candidate hydrocarbon degradaters (Lloyd et al., 2006; Bombach et al., 2009; Abu Laban et al., 2015; Kleindienst and Knittel, 2020). Desulfocapsa was specific to MB.R, while the others were prevalent in CB and MB. The distribution of SRB showed a clear species-specific pattern in correlation with geochemical variables, indicating the diversity of SRB taxonomy and its correlation with the environment (Figure 4B). Sulfurovum was manually divided into three clusters based on the phylogenetic relationships (Supplementary Figure 4). Clusters 1 and 3 were largely specific to MB.R, while Cluster 2 was widely distributed from CB to MV. Clusters 1 and 3 were negatively correlated with nitrate (Spearman, R < −0.3, p < 0.05). Cluster 2 negatively correlated with sulfate and calcium (Spearman, R < −0.3, p < 0.05) but positively correlated with ammonium, silicate, and phosphate (Spearman, R > 0.3, p < 0.05). Milano-WF1B-44 species were prevalent in MV and also showed a negative correlation with sulfate and calcium (Spearman, R < −0.3, p < 0.05) but a positive correlation with ammonium, silicate, and phosphate (Spearman, R > 0.3, p < 0.05, Supplementary Figure 5).

To explore potential carbon and sulfur metabolic pathways and their changes in different aggregates, metagenomic approaches were applied to profile microbial functional genes (Supplementary Table 1). The pmoA/B genes encoding methane monooxygenase and the mcrA gene encoding methyl-coenzyme M reductase were prevalent in MB.R. and MV, but they were specific to bacteria or archaea, respectively, indicating the difference in methane metabolism pathways in these two kingdoms. The aromatic hydrocarbon biodegradation-related genes, bssA and dch (Boll et al., 2002; Winderl et al., 2007; Porter and Young, 2014), were prevalent in CB and MB.R and dominantly affiliated with SRB, indicating the deep involvement of SRB in aromatic hydrocarbon oxidation. The chain hydrocarbon oxidation-related fadA gene, encoding acyl-CoA dehydrogenase in β-oxidation pathway (Dong et al., 2020), was carried by diverse bacterial lineages and also prevalent in CB and MB.R. For carbon fixation, the cbbL gene encoding ribulose-bisphosphate carboxylase for the Calvin-Benson-Bassham (CBB) cycle and the aclB gene encoding ATP-citrate lyase of the reductive tricarboxylic acid (rTCA) cycle were both enriched in MB.R. The cbbL gene was carried out by some unclassified Gammaproteobacteria, implying the Milano-WF1B-44 clade. However, among the aclB gene carriers, we found Campylobacterota and archaea lineages. For sulfur metabolisms, genes encoding the essential enzymes in the dissimilatory sulfate reduction (DSR) pathway, sat and dsrA, were rare in NS but increased in abundance since CB, indicating an increase in DSR metabolism and peaked in MB.R. Other than the known SRB, lineages from Campylobacterota and Gammaproteobacteria were also supposed to carry these genes. The sqr and fccB genes encoding sulfide (quinone oxidoreductase and sulfide dehydrogenase for sulfide oxidation) exhibited a gradual increase from NS to MV. Campylobacterota and unclassified Gammaproteobacteria were carrying these genes, suggesting the potential involvement of Sulfurovum and Milano-WF1B-44 in sulfide oxidation in the cold seep. A similar pattern was also found in the SOX system-related genes, which were dominantly affiliated with SRB, unclassified Gammaproteobacteria, and other bacterial lineages (Figure 5).

We constructed the potential carbon and sulfur metabolic pathways for each aggregate based on the annotated genes. Methane metabolic pathways, including aerobic and anaerobic oxidation of methane (AeOM and AOM), were identified in MB.R and MV. Complete benzoyl-CoA and β-oxidation pathways for aromatic and chain hydrocarbon degradation have been identified since CB, indicating that hydrocarbon degradation potentially could vary in intensity between different aggregates. The CBB cycle was only present in CB and MB for the carbon fixation mechanisms, while the rTCA cycle was specific to MB and MV, showing a noticeable preference for carbon fixation patterns. For the sulfur cycle, DSR, the SOX system, and the sulfide oxidation enzymes were all present since CB, indicating complete sulfur cycling in all aggregates but probably varying in metabolic intensity. None of the pathways mentioned above was identified in NS (Figure 6).