Deep-Sea Carbonates Are a Reservoir of Fossil Microbes Previously Inhabiting Cold Seeps

Carbonates are globally distributed particularly around deep-sea cold seeps. The embedded microbes are fossil records of the past bioprocess but metagenomes of the carbonates have not been fully studied. In this study, we report microbial community structures and genomes of dominant species in cold-seep carbonates from the South China Sea (SCS) and Gulf of Mexico (GoM). The carbonates contained both anaerobic microbes represented by methane oxidizing archaea (ANME) and aerobic ammonia-oxidizing archaea (AOA). The samples from GoM were mostly composed of small microbial groups, indicating heavy degradation of the fossil microbes. The composition of the carbonate communities differed from that of cold seep sediments, suggesting alteration of cold-seep microbial structures during formation and weathering of carbonates. Extraction of 18S rRNA genes from metagenomic reads revealed prevalence of fungal species in the carbonates of the GoM. Genome binning resulted in 10 genomes for dominant prokaryotic species. The ANME genomes showed a short genetic distance to the relatives from the current cold seep sediments; the AOA genomes were affiliated with alpha ecotype dominating deep-sea sediments. Our study reports the genomes in ancient carbonates and sheds lights on microbial role in formation and bioweathering of carbonates.


INTRODUCTION
Along continental coastal slopes, there are cold seeps on the seafloor (Boetius and Wenzhofer, 2013). Carbonate is formed in the deep-sea cold seeps to trap CO 2 into minerals, which is one of important carbon capture processes on Earth. In cold seeps, anaerobic methane oxidizing archaea (ANME) consume methane to generate CO 2 and then secret HCO 3 − (Reitner et al., 2005;Timmers et al., 2017), which is then precipitated with calcium and other metal ions. Syntrophic partners of ANME-1, ANME-2, and ANME-3 are sulfur reducing bacteria (SRB) that provide electron acceptors for the anaerobic oxidation of methane (AOM) (Niemann et al., 2006;Skarke et al., 2014;Li et al., 2020). An archaeal family Candidatus Methanoperedenaceae (known as ANME-2d) may reduce nitrate, nitrite and metal oxides for independent AOM (Niemann et al., 2006). As carbonate deposits are being formed, solution, particles and microbes in the cold seep sediments were embedded into the minerals. Therefore, carbonates can probably serve as an important reservoir of the past geological settings and fossil microbes. In the Black Sea, ANME and SRB were firstly revealed by microscopy (Reitner et al., 2005). Highly diverse microbial communities of deep-sea cold seep carbonates from mud volcano in Mediterranean Sea and cold seep in the Okinawa Though were then revealed by 16S rRNA gene cloning and GeoChip (Heijs et al., 2006;Li et al., 2018). ANME plus other archaeal groups unrelated to methane metabolism were also identified in the carbonates (Heijs et al., 2006). Aside from the anaerobes, there are aerobic microbial groups such as methanotrophic bacteria (Birgel et al., 2011), suggesting that the conditions and positions for the carbonate formation are variable. The precipitation of carbonate may therefore occur in both oxic and anoxic layers of a cold seep as a result of microbial consumption of methane occurring at different oxygen levels (Cavagna et al., 1999). On the other hand, carbonate deposit might be dissolved by sulfur oxidizing bacteria (SOB) (Leprich et al., 2021). This will release the previously buried microbes into the pore water of the cold seep and marine bottom waters, providing additional nutrient source and biodiversity for cold seep ecosystem. Discharging CO 2 into deep-sea water with the dissolvement of carbonates can probably affect the global climate ultimately. However, the SOB have not been detected in carbonate deposits yet. At present, using metagenomic methods for microbial study on cold seep carbonates is still sparse. The genomes of fossil microbes and the community structures in cold seep carbonates at different ages have not been investigated. Whether the fossil ANME and SRB genomes resembled those of the relatives in the current cold seep remains a question.
In the South China Sea (SCS) and Gulf of Mexico (GoM), cold seeps were discovered at numerous sites (Roberts et al., 2010;Feng et al., 2018). Microbial communities and metagenomes have been studied to demonstrate the role of microbes in element cycling and subsequent formation of cold seep ecosystem (Skarke et al., 2014;Li et al., 2020). Biomarkers derived from ANME and bacteria have been examined and mineral components were used to indicate methane flux rate of cold seeps in the past (Bian et al., 2013;Feng and Chen, 2015;Ge et al., 2015). At present, the microbial communities in the cold seeps of the two areas have not been reported and there are not metagenomics data for better elucidation of microbial functions. In this study, we selected four sites in the SCS and GoM to compare their fossil microbes embedded in the carbonates in terms of community structure, diversity and phylogenomic divergence. The age of these carbonate samples had been estimated to be 1.2-46 ka (Birgel et al., 2011;Bian et al., 2013;Feng and Chen, 2015), which allowed to evaluate the impact of lasting time on the preservation of fossil microbes. Our results revealed the presence of a wide spectrum of microbes in the carbonates. The genomes of the dominant species showed a short genetic distance with their current relatives, which further highlights microbial role in carbonate precipitation process and age-dependent deterioration of fossil microbes.

Sampling
During R/V "Dayangyihao" in June 2013, Jiaolong submersible collected carbonate samples from site F of the northern SCS at depth of 1,120-1,143 m. The carbonate sample collection from the SCS was also described elsewhere (Feng and Chen, 2015) (Supplementary Table 1). Carbonate samples were also taken from the northern GoM (Figure 1). GC140 from the Green Canyon (hereafter MEX1) was taken by remotely operated vehicle (ROV) Johnson-Sea-Link dive 2591 at the Green Canyon at depth of 260 m in 1989 (Bian et al., 2013). AC645 (hereafter MEX2) was obtained from the Alaminos Canyon by Alvin dive4197 at depth of 2,200 m (Birgel et al., 2011). The carbonate samples were cut into smaller ones (Figure 1). About 10 g inner parts of the carbonates were digested with 2.5 U DNaseI in 1 L buffered ddH 2 O for 30 min in 25 • C. The carbonates were then cleaned with distilled water under ultrasonication for three times, each lasting for 30 min, followed by an ultrasonication treatment in 70% ethanol for 30 min. The clean carbonates were dissolved in 1% acetate solution. Triplicates were performed for each of the samples.

DNA Extraction, 16S rRNA Gene Amplification and Illumina Sequencing
The carbonate clays after dissolution were used for DNA extraction using a PowerMax Soil DNA isolation kit (Mo Bio Laboratories, Carlsbad, CA, United States) following the manufacturer's instruction. The 16S rRNA gene amplicons of the samples were obtained by a PCR reaction using universal primers 341F (5 -CCTAYGGGRBGCASCAG-3 ) and 802R (5 -TACNVGGGTATCTAATCC-3 ) that target V3-V4 region of 16S rRNA genes (Wang and Qian, 2009). A set of 8-nucleotide barcodes was added to the primers to distinguish the samples (Supplementary Table 2). The PCR was executed as below: denaturation firstly at 98 • C for 10 s; 30 cycles including denaturation at 98 • C for 10 s, annealing at 50 • C for 15 s, and extension at 72 • C for 30 s; and with a final extension at 72 • C for 5 min. The 16S rRNA gene amplicons were prepared for Illumina sequencing by TruSeq R Nano DNA LT Kit (Illumina, San Diego, CA, United States). The amplicon libraries were sequenced with Illumina MiSeq platform (2 × 300 bp). Raw reads were separated for different samples using the barcode information and were then filtered by NGS QC Toolkit (v2.3.3) (Patel and Jain, 2012) with default parameters (base quality <20 for >70 per cent of a given read). The paired reads were merged with PEAR (v0.9.5) . The resulted fragments were processed with QIIME2 pipeline (Caporaso et al., 2010b). The fragments with similarity >97% were firstly clustered into OTUs by UCLUST (Edgar, 2010). The longest fragment in each OTU was selected as the representative read and then used for alignment with the references in the Silva database (version 132) by PyNAST (Caporaso et al., 2010a). Chimeric reads were removed based on the result of ChimeraSlayer (Haas et al., 2011). Taxa assigned to mitochondria, chloroplast, and eukaryotes were deleted. Taxonomy assignment of representative reads was finally conducted by online silvangs 1 . Correlation analysis of genera was carried out by SparCC method (100 permutation).
The 16S and 18S rRNA genes in metagenomics reads (16S and 18S miTags) were predicted using rRNA_HMM (Huang et al., 2009). The reads were assigned to OTUs at 97% similarity level using QIIME v1.9.1 (Kuczynski et al., 2011). Taxonomic classification of the representative 16S miTags in the OTUs was performed using vsearch (Rognes et al., 2016) in reference to SILVA132 (Quast et al., 2013) database. For eukaryotic classification with the same pipeline, those 18S miTags assigned to Metazoa were not taken into account.

RESULTS AND DISCUSSION
The carbonate samples from four sites of the north SCS and GoM ( Figure 1A and Supplementary Table 1) were cementlike minerals composed of nodules, bivalve shells and tiny pores (Figures 1B,C). Previous studies have provided mineral components, age estimate and isotopes of the samples (Birgel et al., 2011;Bian et al., 2013;Feng and Chen, 2015). For the SCS samples, their age ranged between 6.4 and 11 ka; for the GoM samples, the age is largely higher than 10 ka (Birgel et al., 2011;Bian et al., 2013;Feng and Chen, 2015). Triplicates of DNA extraction were executed for different samples from each site. DNA was successfully extracted for 16S rRNA gene amplification and metagenomics work for 10 of 12 carbonate triplicates.
Using Illumina sequencing, we obtained 66,695 16S rRNA amplicon reads for ten samples with a minimum of 2,627 reads (Supplementary Table 2). All the rarefaction curves of the observed OTUs reach a plateau (Supplementary Figure 1), which indicates that the sequencing depth is sufficient to cover the whole microbial communities. A sample from Green Canyon of GoM (MEX1-R1) is associated with the highest number of OTUs (n = 203) and Shannon index of 7.47 (Supplementary Table 2). The carbonates were dominated by Firmicutes, Proteobacteria, and Euryarchaeota (Supplementary Figure 2). However, the microbial communities differed among the samples to some extent. All the carbonate samples contained ANME species such as ANME-1a and ANME-2a-c, which were dominant Archaea in the samples (SCS1: 14.37%; SCS2: 7.39%; and MEX2: 4.87%) except MEX1 (Figure 2). This suggests that the carbonates were formed during methane seeping of the sites. As syntrophic partners of ANME, SRB in comparative percentage were only identified in SCS1-R1 and SCS1-R2. Moreover, ANME-2d that does not rely on SRB was not present in our samples. It is therefore likely that other carbonate samples (such as SCS2) that lacked SRB were formed above sulfur-methane transition zone by other microbes or abiological process (Boetius and Wenzhofer, 2013). The ANME archaea in the carbonates have likely been transferred to near-surface oxic sediment by seeping solution. This was indicated in 6-8 cm layer (JLS_3) of Jiaolong cold seep, where abundant ANME1 archaea were not associated with SRB in comparative amount (Wu et al., 2018). Paucity of syntrophic SRB was also observed in most of in ANMEenriched carbonates formed along coastal slopes of the North America (Leprich et al., 2021). Interestingly, the high relative abundance of ANME archaea in SCS samples was not associated with a comparatively high content of carbonate content in the samples. The previous geochemical analyses showed a carbonate content ranging from 37 to 67% for the SCS samples, much lower than that of 93-96% for the MEX2 sample (Birgel et al., 2011;Feng and Chen, 2015). Therefore, other environmental factors and microbial players might impact on the carbonate content in cold seep.
Ammonia-oxidizing archaea (AOA) represented by Nitrosopumilaceae were prevalent in MEX1, while they were hardly detected in the SCS1 and SCS2 samples. Since AOA were present in oxic subsurface sediment layers and marine waters , carbonate vents were probably formed in anoxic layer and later extruded into shallower oxic sediment layers in the MEX sampling sites. In addition, respiration of heterotrophic microbes in the oxic layer can also generate CO 2 for formation of AOAcontaining carbonates. For the replicates of SCS1, the third sample SCS1-R3 was largely occupied by Bacillaceae and Sphingomonadaceae, suggesting discrepancy between the samples from the same site in the composition of microbial communities. Bacillaceae represented by Bacillus was also ranked as the first and second most abundant genera in SCS2-R1 and SCS2-R2 samples, respectively. In a previous report, several Bacillus strains were experimentally confirmed to be able to dissolve carbonates (Subrahmanyam et al., 2012). Overall, the samples from the GoM were more diversified than those from the SCS as the minor groups occupying <2% of the community were totaled to be more than half of the communities in the GoM samples (Figure 2). The loss of dominant microbial groups in the MEX2 suggests natural degradation of the cemented biomass under the long-term preservation. This sample had been estimated to be formed about 10 ka ago (Birgel et al., 2011). We cannot preclude the possibility that some microbes in the communities might dissolve the carbonate to take the fossil microbes as nutrient. In this study, we have removed the surface carbonate, and therefore most of carbonate-dissolving microbes could not be examined herein.
Correlation analysis revealed relationships between the major microbial groups in the SCS and MEX carbonates (Supplementary Figure 3). High betweenness centrality was found for Nitrosococcaceae and Nitrospinaceae which are known ammonia-oxidizing bacteria (AOB) and nitrite-oxidizing bacteria (NOB), respectively. This seems to be evidence of a microbial network centralized on aerobic nitrogen cycle in the carbonates. Nitrosococcaceae was also found to be a central node of the network between microbial habitants in Karst caves (Zhu et al., 2019), suggesting an unknown role of nitrifying microbes in carbonate formation. The nodes for the AOB and NOB were also negatively related to the ANME species, which agrees with their prevalence in oxic and anoxic conditions, respectively.
Next, we compared the microbial communities of cold seep sediments (JLS) and the carbonates of site F, along with those of MEX1 and MEX2, using principal coordinate analysis (Figure 3). F1 and F2 explain 53% of the differences between the samples. The cold seep sediments were obtained from an active seeping vent. From the surface JLS_1 to bottom JLS_7, four layers of the 14-cm core were studied in our previous work (Wu et al., 2018), demonstrating a transition of SRB to SOB at about 8 cm below surface. In this study, the cold seep sediments were obviously separated from the carbonates, although all were from site F of SCS. This suggests that carbonates cemented a slightly different group of microbes, compared with the cold seep sediments, or the communities have been shifted during carbonate weathering. FIGURE 3 | Principle coordinate analysis of microbial communities at genus level. The Bray-Curtis dissimilarity was calculated using community structures at genus level and was then used for the principle coordinate analysis.
FIGURE 4 | Eukaryotic community in MEX2 using 18S miTags. The 18S rRNA miTags were extracted from MEX2 metagenome. The classification was conducted using QIIME2 against SILVE 132 database for 18S rRNA genes.

Metagenome-Centric Analysis
We obtained a total of 49.2 Gbp metagenome data and 328 million Illumina reads for carbonate samples from the four sites. After quality control, 36.6 Gbp metagenomic data were retained for subsequent work. A total of 104,252 16S miTags were retrieved from the metagenomic reads for detection of prokaryotic community. The primer-independent 16S miTags allow for unbiased detection of microbial communities. Proteobacteria, Firmicutes and Euryarchaeota were still abundantly present in the metagenomes (Supplementary Figure 3). In addition, there are a considerable percentage of the 16S miTags were derived from candidate phyla such as GN04 and Hyd24_12, indicating a huge reservoir of rare prokaryotic species preserved in the carbonates. These rare species had been identified in global cold seep areas (Ruff et al., 2015). We extracted 4,144 18S miTags from MEX2 metagenome, while less than 100 in other metagenomes. A total of 132 18S OTUs were obtained using 97% similarity. After filtration of OTUs for animals, we selected 456 18S miTags for in a brief survey of eukaryotic communities in MEX2. Fungus is the most prevalent group, accounting for 37% of the miTags, followed by SAR superphylum (Figure 4). Most of the fungi were affiliated with Malassezia, Cladosporium, and Aspergillus. About 52% of the 18S miTags were unclassified, suggesting a large number of "dark matter" eukaryotes in deep-sea cold seep. Fungi are the dominant eukaryotic groups in carbonates in Lost City hydrothermal field (López-García et al., 2007). Fungi dominated deep-sea sediments globally and therefore they might be embedded during the carbonate formation (Wang et al., 2014;Rédou et al., 2015). Nevertheless, the presence of abundant fungi in the carbonates also implies the role of fungi in utilization of the fossil microbes as nutrient after the formation of the carbonate as has been argued previously (Bindschedler et al., 2016). Interestingly, the MEX2 samples also harbored a high percentage of minor prokaryotic groups, suggesting degradation of the embedded prokaryotic fossil microbes by fungi in the carbonates. However, this hypothesis warrants future efforts.
Twenty-five MAGs were obtained from the four metagenomes, and most of them (85%) were binned from SCS1 and MEX2. Five of them were contaminants and had been removed. The 20 clean MAGs with completeness higher than 50% were classified with GTDB-tk (Chaumeil et al., 2019; Table 1). The taxonomic novelty could be indicated by RED value that may serve as an estimator of genome novelty in GTDB (Parks et al., 2018). The average of the RED values is 0.80 ± 0.11, which means that most of the MAGs fall into the range of RED values for new families (Parks et al., 2018). The lowest RED value of 0.48 in our MAGs indicates identification of novel class in Gammaproteobacteria.
The MAGs were used to build a phylogenomics tree with the reference genomes of the closest relatives to examine the classification result and to show genetic distance. These MAGs were affiliated with eight phyla (Figure 5). Four of them were assigned to ANME-2c and AOA, which is consistent with our 16S rRNA gene amplicon result. For the ANME MAGs from the site F, they highly resemble those obtained from the cold seep sediment of the same sampling site, regardless of a short genetic distance. In the ANME-2 MAGs, most of the functional genes involved in AOM were identified, including those coding for methyl coenzyme M reductase system (Mcr), anaerobic carbonmonoxide dehydrogenase, CODH/ACS complex (Cdh), 5,10-methylenetetrahydromethanopterin reductase (Mer), tetrahydromethanopterin S-methyltransferase (Mtr) and heterodisulfide reductase (Hdr) (Supplementary Table 3). The products of the AOM process are probably formate, CO 2 and acetate as indicated by the presence of the related genes in the MAGs. The AOA MAGs of MEX1 sample were grouped with Nitrosopumilus maritimus (Figure 5), which belongs to alpha type AOA that was isolated from marine sediment . The biochemical catalysis between CO 2 and HCO 3 − is crucial for the formation of carbonates. ANME-1 and ANME-2 archaeal genomes harbor diversified AOM processes that convert methane into CO 2 , acetate and probably propanoate (Timmers et al., 2017;Li et al., 2020). Carbonate dehydratase is a bidirectional enzyme that acts on a balance between CO 2 and HCO 3 − , which is critical for local pH condition and carbonate stabilization. A total of 83 carbonate dehydratase genes were identified in the assembled contigs of the metagenomes. Five of them were associated with >1× coverage (Supplementary Table 4). With respect to the highest similarity to known homologs, we predicted that these carbonate dehydratase genes were probably derived from ANME-2 (84% identity), Pseudoalteromonas (100% identity), Novosphingobium (99.5% identity), and candidate division FCPU426 (69.9% identity). The finding of the species that encode the carbonate dehydratases likely revealed the major microbial groups that played an important role in formation and/or dissolution of carbonates.
In this study, we detected microbial communities and their representative genomes from carbonate samples of two cold seep areas. Although the two sites were far from each other, they all contained similar microbial consortia of cold seep regardless of different percentages of dominant species.
The difference in the community structure was probably resulted from the discrepancy in carbonate age, cold seep type, and microenvironment. Aside from ANME, there were perhaps unnoticed microbial inhabitants mediating carbonate formation in oxic layer above methane-sulfate transition zone. Dissolution and reprecipitation of carbonate under effect of both biological and chemical processes on an existing carbonate might also have changed the embedded microbial communities. More samples with clear geochemical background may help explain the variation of fossil communities in carbonates and highlight their potential significance to global climate change.

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

AUTHOR CONTRIBUTIONS
YW conceived and designed the study and wrote the manuscript. WL performed the experiment and analysis. QL, YZ, and ZG provided the analysis pipelines. DF provided the samples and associated metadata. All authors agree with the final version of the manuscript.