The sediment sample core-10 was collected by means of push core at a deep-sea vent site at the Guaymas Basin in the Gulf of California in October 2008, during Alvin dive 4460, which was part of cruise AT15-25. The sediment used in this study was oily and covered with a white microbial mat. The temperature, as measured by the Alvin temperature probe, increased from 7 to 10 °C at the surface to 31 °C at 6 cm below the surface (cmbsf). Two sediment layers, at 4–6 and 8–11 cmbsf, were designated as M8 and M10, respectively, and used for DNA isolation. The methods used for DNA extraction, metagenome library construction, and sequencing are described in Ref. [4].

The raw reads from metagenomic sequencing were trimmed with Sickle v1.33 (https://github.com/najoshi/sickle) and assembled using IDBA-UD v1.1.1 [15] with the following parameters: “mink 52, maxk 92, step 8.” The coverage was determined by recruiting raw reads to assembled contigs using Bowtie v2.2.8 [16] with the parameter “–very-sensitive,” and was calculated using SAMtools v1.3.1 [17]. Binning of assembled sequences longer than 5 kb was performed using tetra-nucleotide frequencies, guanine-cytosine (GC) content, and sequencing coverages using MetaBAT v0.32.4 with the parameter “–specific” [18]. The retrieved MAGs were further checked manually using ESOM v1.1 [19] and mmgenome package v0.6.3 [20]. The completeness and contamination of each MAG were estimated using CheckM v1.0.7 with the “lineage_wf” function [21].

In addition to the MAGs retrieved from the present study, 44 reference bathyarchaeotal genomes were downloaded and used in the current analysis. These MAGs were retrieved from samples from Aarhus Bay sediment [2], Surat Basin coal-bed methane wells [12], Colorado soil samples [22], Guaymas Basin hydrothermal sediments [4,23], Northern California grassland [24], White Oak River estuarine sediments [9], Juan de Fuca basalt-hosted fluids [25], and other globally distributed sample sites [26].

This Whole Genome Shotgun project has been deposited at GenBank under the accession number PRJNA418890.

A total of 3738 reference bathyarchaeotal sequences longer than 900 nt were collected from the SILVA 132 database [27] and clustered into 1547 operational taxonomic units (OTUs) using CD-HIT [28] with the parameter “–c 0.97-G 1.” Another 23 16S rRNA gene sequences were predicted from the MAGs mentioned above using CheckM [21]. All of these 16S rRNA gene sequences, together with reference sequences from Thaumarchaeota and Aigarchaeota which were used as outgroups, were aligned using MAFFT v7.313 [29] with the parameters “–localpair” and “–maxiterate 100.” The phylogenetic tree was generated using maximum likelihood analysis in RAxML v8.2.8 [30] with the GTRGAMMA model, and visualized using Interactive Tree Of Life (iTOL) [31]. Affiliations of bathyarchaeotal sequences to each subgroup were carried out using reference sequences [1,8] as phylogenetic anchors.

Protein coding sequences (CDSs) were determined using Prodigal v2.6.3 [32] with the parameter “–meta.” Functional information for each CDS was collected by comparing amino acid sequences against the National Center for Biotechnology Information (NCBI) “non-redundant” (nr), Kyoto Encyclopedia of Genes and Genomes (KEGG) [33], and Clusters of Orthologous Group (COG) databases [34] with E-values < 10⁻⁵, and was inspected manually. Those predicted CDSs were further searched against the MEROPS [35] and CAZy [36] databases with E-values < 10⁻⁵ to find genes related to the degradation of peptides and carbohydrates. The extracellular peptidases were further confirmed based on the identification of transport signals using PRED-SIGNAL [37].

To identify relevant phylogenetic homologs of bathyarchaeotal MAGs, a total of 60 bathyarchaeotal genomes (16 from the present study and 44 from public databases) were included in the phylogenetic analysis, along with other archaeal reference genomes. Following a previous study [9], 16 conserved ribosomal genes (ribosomal protein L14, L15, L16/L10AE, L18, L2, L22, L24, L3, L4, L5, L6P/L9E, S10, S17, S19, S3, and S8) were considered to undergo no horizontal gene transferring [38], and were used to infer the phylogenetic affiliation of the inspected bathyarchaeotal MAGs. Homologous sequences were individually aligned using MAFFT v7.313 [29] with the parameters “–localpair” and “–maxiterate 100,” and then concatenated in series. The phylogenetic trees were generated using maximum likelihood analysis in RAxML v8.2.8 [30] with the PROTGAMMA model, and visualized using iTOL [31].

Ancestral genome reconstruction was implemented in the Count software [39]. The phyletic patterns (gene family presence/absence and gene copy number) of hydrothermal-adaptation-related genes were mapped to the bathyarchaeotal phylogeny based on the 16 concatenated conserved ribosomal protein sequences mentioned above. The phyletic pattern was visualized using iTOL [31] and plotted manually.

In total, 16 bathyarchaeotal MAGs were recovered from the Guaymas Basin M8 and M10 metagenomic datasets (Table 1). Together with 44 published bathyarchaeotal genomes, 23 16S rRNA gene sequences were predicted (Appendix A. Table S1) using CheckM [21]. Furthermore, an additional 1547 representative bathyarchaeotal 16S rRNA gene sequences with a maximum 97% identity cutoff were downloaded from the SILVA databases and used as references. A 16S rRNA gene phylogenetic tree was constructed (Fig. 1(a) and Appendix A. Fig. S1), whose topology was highly similar to and coherent with results from previous studies [1,4,8,14,40]. To be specific, Bathy-21 and -22, which were located at the base of the tree, branched monophyletically from the other subgroups. These two subgroups, together with subgroup Bathy-23, were recently named by Zhou et al. [8]. The related 16S rRNA gene sequences are listed in Appendix A. Table S2. The phylogeny for most bathyarchaeotal subgroups are relatively stable, with variations in only a few branches’ orders still remaining due to the use of different reference sequences, methods, and parameters in different studies [1,4,8,14,40]. Instability in the tree topologies is probably due to insufficient sequences for phylogenetic analyses, indicating the extremely high diversity of the phylum Bathyarchaeota.

A phylogenomic tree of the MAGs from the current study (16) and public database (44) was further constructed using 16 concatenated conserved ribosomal protein sequences (see Section 2.5 and Fig. 1(b)). A total of 54 MAGs were classified into 11 different subgroups (Table S1), whereas six MAGs (ex4484_135, ex4484_205, M10_bin214, M10_bin241, M8_bin163, and JdFR-01) were categorized as unclassified due to the lack of 16S rRNA gene sequences. Similar to the phylogenetic analysis based on 16S rRNA gene sequences, the subgroups Bathy-21 and Bathy-22 are still located at the root of the phylogenomic tree, which suggests that these two linages may represent the ancient Bathyarchaeota types. Subgroups 1, 3, 5, 6, 8, 15, and 17 have corresponding genomes and show similar topology with the 16S rRNA gene phylogenetic tree. Nevertheless, the phylogeny for subgroups 13 and 23 was not matched between the 16S rRNA-based and concatenated ribosomal protein-based phylogenomic analyses (Table S1). Considering the higher resolution and stability of the phylogenetic analysis based on concatenated conserved ribosomal protein sequences compared with 16S rRNA genes, the phylogenomic tree was used for further analysis.

Among 60 bathyarchaeotal MAGs, six MAGs (ex4484_135, ex4484_205, M10_bin214, M10_bin241, M8_bin163, and JdFR-01) remained in an unclassified subgroup, while three MAGs (M10_bin1903, DG-45, and RBG_13_46_16b) contained potentially high contamination > 10%; these nine MAGs were not used for the core and specific metabolic analysis of different subgroups. The remaining 51 phylogenetically affiliated MAGs formed 10 bathyarchaeotal subgroups, and their metabolic potentials were systematically compared (Fig. 2 and Table 2). While Bathy-3 had only one representative MAG, other subgroups contained two or more MAGs. Shared core metabolic features among different lineages were discovered, including the pathways for the utilization of detrital proteins, lipids, and benzoates. Detrital proteins are considered to be one of the most abundant components of marine organic matter [41]. Pathways involved in the degradation of proteins were present in all bathyarchaeotal subgroups (Appendix A. Table S3). The detected extracellular peptidases mainly belonged to the serine, cysteine, and metallo peptidase families (Table S3), while subtilisin (S08A), papain (C01A), and gamma-glutamyl hydrolase (C26) existed in most bathyarchaeotal subgroups. These observations further agreed with previous studies [2,4], suggesting that detrital proteins and/or amino acids are used as major carbon and energy sources for Bathyarchaeota. The degradation of lipids and benzoates was another core genomic feature. Except for one subgroup (Bathy-21 for lipids and Bathy-3 for benzoates), the majority of the bathyarchaeotal subgroups contained long-chain acyl-CoA synthetase (fadD) and acylphosphatase (acyP, Table 2 and Appendix A. Table S4), the key genes for lipid and benzoate degradation, respectively. In general, these metabolic potentials indicate the capability of Bathyarchaeota to degrade diverse complex organic compounds for an organoautotrophic life strategy, as was recently suggested [11].

The processes for energy conservation were mostly identified as the Embden–Meyerhof–Parnas (EMP) glycolysis pathway and the reductive acetyl-CoA (Wood–Ljungdahl, WL) pathway (Fig. 2 and Table 2). Most of the Bathyarchaeota genomes contain near-complete enzymes for the utilization of glucose to produce acetyl-CoA by the glycolysis process (Table 2 and Table S4). Although the pyruvate kinase gene (pyk), which converts phosphoenolpyruvate (PEP) into pyruvate in the last step of EMP, was found in only six out of 10 subgroups, another gene encoding phosphoenolpyruvate carboxylase (ppc), which converts PEP into oxaloacetate, was present and could be used for an alternative pathway.

Bathyarchaeota can also use the WL pathway to generate energy. In nearly all subgroups, most of the carbon metabolic pathways lead to acetyl-CoA production—including the autotrophic WL pathway which uses CO₂ as the carbon source. All the analyzed subgroups possessed an archaeal-type WL pathway for autotrophic inorganic carbon fixation using tetrahydromethanopterin (H₄MPT) as the C₁-carrier, which is prevalently found in Archaea (Table 2 and Table S4). Interestingly, a bacterial type of the WL pathway that uses tetrahydrofolate as an intermediate was also found in nearly all genomes from the Bathy-1, -15, and -17 subgroups. Although gene encoding formate dehydrogenase (fdh), which converts CO₂ to formate, was not found in these three subgroups, all other genes essential for the bacterial WL pathway were present, including carbon-monoxide dehydrogenase (coo), formate–tetrahydrofolate ligase (fhs), methylenetetrahydrofolate dehydrogenase (fol), methylenetetrahydrofolate reductase (met), and acetyl-CoA decarbonylase (cdh). Therefore, these three subgroups might contain unidentified enzymes for the conversion of CO₂, using both H₄MPT and tetrahydrofolate as the C₁-carrier.

Metabolic predictions based on the genomic information of different subgroups inferred group-specific metabolic properties that are not commonly shared in all lineages (Fig. 2, Table 2, and Table S4); these pathways include carbohydrates utilization, the tricarboxylic acid (TCA) cycle, acetogenesis, and sulfur-related metabolism. Different preferences of bathyarchaeotal subgroups for uptaking and degrading a wide variety of carbohydrates were found (Table 2 and Table S4), including cellulose, glycogen, galactan, and mannan. A complete TCA cycle was only detected in the subgroups Bathy-1, -15, and -17; other subgroups such as Bathy-6, -8, and -22 possessed incomplete TCA cycles, probably due to the incompleteness of the genomes. All 10 MAGs belonging to Bathy-23 (as well as the highly contaminated MAG M10_bin1903) lacked the citrate synthase (cs) and aconitate hydratase (acn) used in the first and second step of the TCA cycle, as well as succinyl-CoA synthetase (suc), which converts succinyl-CoA into succinate. The results above suggest that these widespread Archaea might be heterotrophs with extensive abilities to efficiently use organic matter—a survival strategy in carbon-starved deep marine sediments. This suggestion is consistent with the results of previous stable isotopic probing experiments on ¹³C-organic compounds, as well as with the correlation between bathyarchaeotal abundance and the content of total organic carbon (TOC) in diverse habitats [5,42–44].

Acetogenesis is an important feature shared by some bathyarchaeotal subgroups (Table 2 and Table S4): a biological process in organic carbon cycling in subsurface sediments [45]. Subgroups Bathy-15 and -17 have the potential to carry out phosphorylation catalyzed by adenosine diphosphate (ADP) forming acetyl-CoA synthase (acd), while Bathy-1, -21, -22, and -23 can also perform acetogenesis via phosphate acetyltransferase (pta) and acetate kinase (ack). The latter two genes were found to be widely distributed in bona-fide bacterial homoacetogens, which produce acetate solely from hydrogen and CO₂. However, further genomic evidence is still required to answer the question of whether Bathy-3, -6, -8, and -13 have acetogenesis ability. Recent findings suggest that acetate is considered to play a crucial role in mediating the microbial carbohydrate utilization, fermentation and respiratory pathways in the Guaymas Basin sediments [23]. In the present study, no less than six out of 10 subgroups of Bathyarchaeota were considered to be potential acetogens, and may play a vital role in carbon cycling in marine sediments.

The genomic evidence also indicated diverse metabolic potentials in different subgroups regarding sulfur metabolic pathways (Table 2 and Table S4). Bathy-1, -21, and -23 contained genes related to the hydrolysis of sulfate esters, which is responsible for the decomposition of organic sulfur compounds. It is notable that these subgroups also contained genes encoding the sulfhydrogenase gene (hyd, also found in Bathy-17 and -21) involved in sulfide formation from polysulfide (sulfidogenesis), and the sulfate transport gene (cysU). These results suggest a formerly ignored role of Bathyarchaeota in global sulfur cycling.

The majority of the bathyarchaeotal subgroups encoded diverse molecular chaperones, such as heat-shock protein (htpX) and small heat-shock proteins (sHSP, Table S4), which were previous considered to express during the heat-shock process [46]. Interestingly, Bathy-21 and -23, which were retrieved from Guaymas Basin hydrothermal sediments, were devoid of the bacterial-type DnaK-DnaJ-GrpE chaperone system that is prevalent in other subgroups (Table S4). This simplified heat-shock system was previously reported to be one of the characteristics of hyperthermophilic Archaea [46], and thus indicates a potential hyperthermophilic feature of these bathyarchaeotal subgroups. Genes encoding reverse gyrase, a marker gene for hyperthermophiles, were detected in Bathy-21, -22, and -23 (Table S4), and in three unclassified bathyarchaeotal MAGs (ex4484_205, ex4484_135, and B242). The ancestral genome analysis (Fig. 1(b)) showed that the putative Bathyarchaeota ancestral lineage possessed the reverse gyrase gene. However, the hyperthermophilic bathyarchaeotal MAGs were exclusively retrieved from Guaymas Basin hydrothermal sediment. Other bathyarchaeotal subgroups may have lost the hyperthermophilic ability during evolutionary adaptation to mesothermal environments. Based on existing bathyarchaeotal MAGs, we hypothesize that the ancient Bathyarchaeota type might be hyperthermophiles, which gradually adapted and evolved into ambient environments. However, more bathyarchaeotal genomic data are needed to further confirm this assumption. The archaeal phylum Bathyarchaeota belongs to the TACK superphylum (or “Proteoarchaeota”). Although many lineages from the TACK superphylum (e.g., Korarchaeota, Crenarchaeota, and Aigarchaeota) have been reported to be hyperthermophilic [47], no bathyarchaeotal MAGs or subgroups have been found to be potentially hyperthermophilic to date. The above results indicated that the subgroups Bathy-21, -22, and -23, which were retrieved from sediments with a high temperature, may be hyperthermophiles. Furthermore, a global survey of the distribution of these subgroups indicated that the majority of these subgroups are found in high-temperature environments (Table S2), which provides geographical evidence of the hyperthermophilic nature of these bathyarchaeotal subgroups.