Introduction

Marine subsurface sediments are enriched with diverse groups of archaeal lineages, although their distributions, ecological roles, and adaptation strategies are understudied [1,2,3,4]. Metagenomic sequencing and single cell genomics have enabled the discovery of novel organisms, the elucidation of new metabolisms, the expansion of known lineages, and the redefinition of portions of the tree of life [5,6,7,8,9,10,11,12,13]. However, there is still a paucity of genomes resolved from the deep marine subsurface, meaning niche specific adaptations in deep biosphere are not yet well understood.

Recently, the tree of life has been expanded with the discovery of the Asgard superphylum, a deeply-branching monophyletic group thought to be some of the closest relatives to the eukaryotic branch of life [10, 14]. Genome analyses of Asgard archaea noted high abundances of eukaryotic homologs and have suggested diverse metabolic functions extending from an autotrophic lifestyle, primarily dependent on carbon fixation via Wood–Ljungdahl pathway and acetogenesis, to a heterorganotrophic lifestyle consuming proteins and aliphatic hydrocarbons, using methyl-CoM reductase-like enzymes [15,16,17]. Asgard members have been proposed to be engaged in symbiotic partnerships involving syntrophic transfers of hydrogen and electrons following the “reverse flow model” [15].

Members of another archaeal phylum in the deep subsurface, Bathyarchaeota, are characterized by their wide metabolic range enabling heterotrophic scavenging of proteins, carbohydrates, short chain lipids, and other reduced compounds as substrates as well as their methane-metabolizing potential [18]. However, bathyarchaeotal genomes have also suggested the potential for carbon fixation and acetogenesis [19]. The evolutionary path describing the acquisition of both methanogenesis and acetogenesis pathways in Bathyarchaeota remains unresolved [19, 20].

Deep sediment from the Costa Rica (CR) margin subseafloor, sampled during the International Ocean Discovery Program (IODP) Expedition 334, was recently shown to host abundant archaea [21]. Here we examine in detail the archaeal genomes recovered from metagenomes from the CR margin and compare them to their relatives recovered from shallower sites. In this study, we report 31 archaeal metagenome-assembled genomes (MAGs) belonging to six different archaeal lineages (Lokiarchaeota, Thorarchaeota, Heimdallarchaeota, Bathyarcheota, Thermoplasmatales, and Hadesarchaea). We thoroughly analyze representative MAGs of two novel archaeal lineages belonging to Lokiarchaeota and Bathyarchaeota. Our analysis suggests that Lokiarchaeota genomes encode for genes to process and degrade aliphatic and aromatic hydrocarbon anaerobically. We also describe a novel Bathyarchaeota genome, which lacks methyl coenzyme M reductase (MCR) complex and possesses a noncanonical acetogenic pathway linking methylotrophy to acetogenesis via the methyl branch of Wood–Ljungdahl pathway. Lastly, we integrate genomic and thermodynamic modeling to underline the ecological and physiological conditions that could drive the syntrophic interactions among CR_Asgard and use phylogenetic inference to elucidate the development of noncanonical acetogenesis in Bathyarchaeota in CR deep sediments.

Materials and methods

Site information and sample collection

Samples used in this study were collected under aseptic conditions from Sites U1378 and U1379 of CR margin during IODP Expedition 334. The sample depths are in the range of 2–93 m below seafloor (mbsf). Detailed site descriptions were previously reported in the IODP Proceedings for Expedition 334 [21, 22].

DNA extraction and sequencing and metagenome-assembled genome binning

DNA extraction and metagenomic sequencing have been described previously [21] and data have been deposited at NCBI GenBank SRA under project PRJEB11766. Metagenomic reads were quality trimmed using Nesoni following default parameters and applying q20 for quality score (www.vicbioinformatics.com/software.nesoni.shtml). Quality-controlled reads in individual samples were assembled separately using IDBA-UD [23] with default settings. Contigs longer than 1 kb were binned into MAGs using MaxBin V2.2.7 [24], further curated manually using VizBin [25], through filtering outlier scaffolds not falling within the same GC%, and differential coverage levels across different CR datasets. The quality and completeness of the MAGs were assessed using CheckM (v.1.0.7) [26]. MAGs from all five metagenomes were dereplicated using dRep (version v2.0.5 with ANI >99%) [27] and the most complete MAG per taxon was selected for downstream analyses. Phylogenetic placements of archaeal MAGs were further determined by Phylosift [28] analysis based on single copy marker genes and phylogenomic analysis using 16 concatenated ribosomal proteins (see description below). Assembled contigs larger than 1 kb were annotated using PROKKA [29]. Encoded proteins were predicted using Prodigal v2.6.3 with the default translation table (Table 11) was applied [30].

Concatenated ribosomal protein phylogeny

A maximum-likelihood tree was calculated based on the concatenation of 16 ribosomal proteins (L2, L3, L4, L5, L6, L14, L15, L16, L18, L22, L24, S3, S8, S10, S17, and S19) using IQ-Tree (v1.6.6) [31] (located on the CIPRES web server) [32]. References sequences used were collected [11] with the addition of more representatives of the Bathyarchaeota and Asgard lineages collected from all MAGs available in JGI_IMG and NCBI database. Evolutionary distances were calculated based on best fit substitution model (VT + F + R10), and single branch location was tested using 1000 ultrafast bootstraps and approximate Bayesian computation [31, 33], branches with bootstrap support >80% were marked by black circles.

Metabolic reconstruction and functional annotation

Predicted proteins from all MAGs were screened using HMMsearch tool against custom HMM databases representing the key genes for specific metabolic pathways [34]. The completion of the pathways was assessed through querying the predicted proteins against KEGG database using BlastKoala tool [35]. Carbohydrate-active enzymes were identified using dbCAN-fam-HMMs (v6) database [36]. Cellular localizations of predicted proteins were identified using SignalP v5.0 [37]. Proteases, peptidases, and peptidase inhibitors were identified using USEARCH-ublast tool [38] against the MEROPS database v12.1 [39]. Transporters were identified using USEARCH-ublast tool [38] against the TCDB database [40]. ESPs were detected using Interpro v75.0 [41]. Phylogenetic distributions of the predicted proteins in each bin were detected through comparing these proteins against the NCBI (nr) protein database via the DarkHorse software [42].

Functional protein-based trees

All functional protein-based trees were built by aligning the query protein sequences to the reference sequences belonging to the same protein family using Muscle v3.8.31 [43]. Reference sequences were collected from AnnoTree [44] using the corresponding KEGG entry as search keyword. Aligned sequences were manually curated using Geneious v9.0.5 (https://www.geneious.com). The phylogenetic trees were computed using IQ-TREE (v1.6.6) [31], through the CIPRES web server [45] and the evolutionary relationships were described using the best fit model. Branch locations were tested using 1000 ultrafast bootstraps and approximate Bayesian computation [31, 33].

Thermodynamic calculations

Initial Gibbs free energy (ΔrG′) calculations were performed for the redox reactions proposed in the CR Lokiarchaeota and their potential partners using eQuilibrator [46] applying pH 8, similar to the approach reported in [47], and reactant concentrations of 1 mM each.

We further confirmed if the coupled redox reactions proposed for the novel CR Lokiarchaeota genomes were feasible, by calculating the Gibbs free energy of five redox reactions (Table S1) under the near in situ conditions of CR sites U1378 and U1379, following the method described in ref. [48]. Gibbs free energy was calculated using the equation:

$$\Delta G_r \,=\, \Delta G_r^0 \,+\, {\mathrm{RT}}\,{\mathrm{In}}\,Q_r,$$

where \(\Delta G_r^0\) and Qr refer to the standard molar Gibbs energy and the reaction quotient of the indicated reaction, respectively, R represents the gas constant, and T denotes temperature in Kelvin. In this study, \(\Delta G_r^0\) was calculated using the thermodynamic data of standard Gibbs free energy of formation of each species and corrected to near in situ pressure and temperature (4 °C), using the R package CHNOSZ [49]. Qr stands for the reaction quotient, which can be calculated with the relation

$$Q_r \,=\, {\Pi}\left( {{\mathrm{a}}_i^{\nu i}} \right),$$

where ai is the activity of species i and νi is its stoichiometric coefficient. ai is the product of chemical species concentration [i] and its activity coefficient γi, which was computed as a function of temperature and ionic strength by using an extended version of the Debye–Huckel equation [50]. Because most of the reactant concentrations were hard to measure or were below detection limits, we assumed 0.1 µM for the concentration of NO3, NO2, SO32−, and H2S. We benzoate (C7H5O2) concentrations in the range of 0.0001–100 µM to explore the feasibility of these reactions over a wide range of substrate concentrations.

Results

Diversity of the dominant archaeal lineages based on 16S rRNA gene analyses

High abundances of archaea across sediment samples collected from five depths of the CR Margin (2, 32, and 93 mbsf from Site 1378 and 22 and 45 mbsf from Site 1379 of IODP Expedition 334), were previously reported [21]. Metagenomic data from these samples were assembled and examined for small subunit ribosomal genes. A total of 126 16S rRNA gene sequences were recovered from the metagenome assemblies. Sequences affiliated to archaea (31 sequences representing 25% of the total) represented a wide diversity of lineages, including Bathyarchaeota (14 sequences representing 11% of the total 16S rRNA gene sequences), Thermoplasma (5 sequences representing 4%) and Lokiarchaeota (4 sequences representing 3%) (Fig. S1).

General genomic features of the abundant archaeal lineages

Across the five depths analyzed, 31 different draft archaeal MAGs were recovered (Table 1). Genomic analyses were only performed on the 11 MAGs showing high completeness and low contamination (above 60% completeness and below 5% contamination). Overall, completeness varied from 32 to 99% with an average of 50% and contamination varied from 0 to 10% with an average of 8% (Table 1). Genome qualities were further assessed by comparing their predicted proteins against the NCBI (nr) database to evaluate the extent of phylogenetic consensus within the binned genomes (Fig. 2a). Overall, the taxonomic affiliations of the majority of the predicted proteins (60–75%) in each MAG agreed with their respective phylogenetic group, except CR_12 in which only 15% of the encoded proteins assigned to Heimdallarchaeota and was excluded from further analysis.

Table 1 Details of the MAGs constructed from the Costa Rica margin.

Phylogenetic placement of the draft genomes was determined using 16 ribosomal proteins (Fig. 1) [11]. The recovered MAGs were affiliated to six different phylogenetic lineages, namely, Lokiarchaeota, Thorarchaeota, Heimdallarchaeota, Bathyarcheota, Thermoplasmatales, and Hadesarchaea.

Fig. 1: Phylogenetic placement of the Costa Rica archaeal draft genomes (genomes of this study are highlighted in red).
figure 1

The maximum-likelihood phylogenetic tree was calculated based on the concatenation of 16 ribosomal proteins (L2, L3, L4, L5, L6, L14, L15, L16, L18, L22, L24, S3, S8, S10, S17, and S19) retrieved from the Costa Rica archaeal genomes and 231 reference archaeal genomes representing 13 different archaeal phyla. The relationships were inferred using the best fit substitution model (VT + F + R10) and nodes with bootstrap support >80% were marked by black circles. Scale bar indicates substitutions per site. The tree is available with full bootstrap values in Newick format in the Supplementary Data.

Distribution of eukaryotic signature protein (ESP) homologs in Costa Rica archaea MAGs

All archaeal MAGs recovered from CR sediments contained eukaryotic signature proteins (ESPs). These ESPs include protein homologs dedicated for information processing, trafficking machinery, ubiquitin system, cell division, and cytoskeleton formation. The Asgard archaea genomes recovered from CR, e.g., Heimdallarchaeota and Lokiarchaeota MAGs (CR_06, CR_11), showed high numbers of eukaryotic homologs and covered broad classes of ESPs (Fig. 2b). However, ESPs were also detected in the Bathyarchaeota, Hadesarchaeota, and Thermoplasmata MAGs recovered from CR (Fig. 2b). Yet, only the Asgard genomes have homologs for cell division and cytoskeleton.

Fig. 2: Phylogenetic distribution and the occurrence of eukaryotic-specific proteins in archaeal MAGs.
figure 2

a Phylogenetic distributions of the predicted proteins encoded by the archaeal MAGs, the phylogenetic assignments were performed through comparing the proteins to the NCBI (nr) proteins database via the DarkHorse software [42]. b Heat map shows the presence (black)/absence (white) patterns of eukaryotic homologs in the recovered archaea MAGs. The color along the left side shows the completeness of the genome bin, as less complete bins would be expected to contain fewer homologs. Only the Asgard archaea contain cell division/cytoskeleton homologs.

Genomic evidence of CR_Asgard hydrocarbon utilization

We screened all of the MAGs, assembled metagenomes, and raw reads for the presence of genes encoding for MCR complex to assess the role of the archaeal members in mediating the degradation of hydrocarbon compounds (Fig. 3a) that are abundant in CR margin [51]. MCR complex genes (mcrABCDG) were completely absent from both the MAGs and the entire metagenomes, which suggests that short chain alkanes are not oxidized using MCR complex in CR sediments. All metagenomic reads and MAGs were screened for possible alternative hydrocarbon degradation pathways using custom HMM searches that specifically targeted key metabolic genes for aliphatic and aromatic hydrocarbon degradation pathways (Fig. 3a). Multiple pathways were successfully identified including glycyl-radical enzymes related genes coupled with n-alkane succinate synthase (AssA) and benzylsuccinate synthase (BssA), which activates n-alkanes and mono-aromatic compounds, respectively, by forming C–C bond between these compounds and fumarate to form hydrocarbon adducts [45, 52]. Benzylsuccinate synthesis is the initial step for aromatic hydrocarbon mineralization, in which benzylsuccinate is converted to benzoyl-CoA [52]. Interestingly, the potential for ATP-dependent Benzoyl-CoA reductase (BCR) complex utilization was identified in some CR Asgard members (Lokiarchaeota, Thorarchaeota, and Heimdallarachaeota) (Fig. 3). This reaction can dearomatize the benzoyl-CoA to dienoyl-CoAs as the first step in aromatic hydrocarbon degradation and then couple this reaction with a beta-oxidation pathway to ultimately produce acetyl-CoA, similar to the mechanism previously reported in the denitrifying bacteria Thauera aromatica [53]. Since Lokiarchaeota CR_06 had the least contamination levels <2%, we used this MAG to verify the presence of aromatic hydrocarbon degradation function mediated by BCR complex in Asgard phyla. Phylogenetic analysis of the BCR subunit B recovered from Lokiarchaeota CR_06 showed their affiliation to class Bzd, which is composed of four subunits (BzdONPQ) (Fig. 3c). This BCR type was originally discovered in Betaproteobacteria, Azoarcus evansii, and the ones detected in CR_06 are closely related to BCRs detected in different β, δ proteobacteria, and other archaeal lineages (e.g., Lokiarchaeota, Bathyarchaeota, and Archaeoglobus) (Fig. 3c) [54]. Subunits P and Q have ATP binding/ATPase functional domains and are linked together by a [4Fe–4S] cluster. Reduced ferrodoxins transfer electrons to the CoA ester-binding domains (subunits O and N), which catalyzes the cleavage of benzoyl-CoA aromatic ring and yields dienoyl CoA product (Fig. 3a) [54]. In CR_06, the four subunits of BCR complex were located in one contiguous operon (CR_06-contig-100_3495), however due to the high fragmentation of the genomes no reliable phylogenetic marker genes were found in the adjoining genomic neighborhoods. A plausible explanation for the high fragmentation levels of the CR_06 MAG (number of scaffolds = 1110), even though the genome exhibited high coverage levels >200×, is that there are high levels of intralineage strain heterogeneity. This is confirmed by the ANI values (ANI = 70–80%) between MAGs of the same phylogenetic group and disparities in their coverage levels (Tables 1 and S3).

Fig. 3: Lokiarchaeotal metabolism and degradation of hydrocarbons.
figure 3

a Metabolic reconstruction of the CR_Lokiarchaeota bin (CR_06). Central metabolic pathways are shown in gray boxes, carbon fixation pathways (WL and rTCA cycles) are shown in pink, electron transport chain (ETC) proteins are shown in yellow, fermentation products are shown in blue boxes, amino acids are shown in brown boxes, enzymes and enzyme complexes are shown in orange circles, energy carriers are shown in red, energy molecules are shown in dark green, metabolite and amino acid transporters are shown in light green, and cation transporters are shown in dark blue. b Distribution patterns of the hydrocarbon degradation pathways among the constructed CR_Archaea genomes. Genomes were clustered based on presence (black)/absence (white) profiles using Euclidean distance and average linkage method. X-axis represents enzymes included in the analysis, bottom bars are the pathways that these enzymes represent, and Y-axis includes bin names/phylogenetic affiliations. c Maximum-likelihood tree of the benzoyl-CoA reductase subunit B. The tree was calculated using the best fit substitution model (VT + F + R7) that describes the evolutionary relationships between BCRB families. The tree was made using reference sequences under the KEGG entry (K04113) collected from AnnoTree [44] and branch location was tested using 1000 ultrafast bootstraps and approximate Bayesian computation, branches with bootstrap support >80% were marked by black circles. Blue and red clades highlight sequences belong to Bzd_O and BCR_B subfamilies, respectively. Scale bar indicates substitutions per site. Sequences from CR_Lokiarchaeota bins was marked with red circles, CR_Thorarchaeota bin was marked with light blue, and CR_Heimdallarchaeota was marked with green. Candidate novel clades present in the Costa Rica metagenomic datasets were marked with purple stars. The tree is available with full bootstrap values in Newick format in the Supplementary Data.

The affiliations of the BCR complex related proteins, BCRA and BCRB in Lokiarchaeota CR_06 were confirmed by the following observations: (1) contigs encoding for the four subunits of BCR complex fall within the same coverage and GC% range of the rest of the contigs in the same MAG (Fig. S2); (2) phylogenetic trees of the BzdO and Q protein sequences detected in the CR_06 were placed as siblings to BCR sequences belonging to Lokiarchaeota MAGs recovered from deep sediment from south China Sea (GenBank accession numbers TET61664 and TKJ21612) (Figs. 3c and S3) [55].

Fate of the degraded aromatic hydrocarbons and syntrophic interactions

CR_ Lokiarchaeota (CR_06) as well as other Asgard archaea in the CR margin could potentially grow heterotrophically on aromatic hydrocarbons. Though, the absence of genes encoding for different types of cytochrome oxidase and anaerobic respiration from the genome content of the Asgard MAGs (CR_06, 07, 08, 11, and 12) indicate their inability to completely mineralize these hydrocarbons to CO2 and H2O, which would allow a high energy yield. They encode for the genes mediating the fermentation of these organic macromolecules to acetate and other reduced products. However, this is thermodynamically unfavorable under CR conditions with a positive ΔrG′ value (ΔrG′ = 196.3 [kJ/mol]) (Fig. 3 and Table S1), suggesting that these genes are selectively maintained. Due to the genome incompleteness (86% complete), we cannot rule out the possibility for the presence of complete aromatic hydrocarbon mineralization pathways using one or more oxidized substrates as electron sinks (Table 1). More likely, however, Lokiarchaetoa are gaining energy through syntrophic interactions with partners capable of oxidizing the biodegradation intermediates.

Here, we inferred the identity of the potential syntrophic partners under marine subsurface conditions by comparing all the possible metabolic and thermodynamic scenarios, gauging each scenario based on the presence of the metabolic pathways in our metagenomic datasets and the thermodynamic feasibility under each condition (Fig. 3b, Tables S1 and S6). We calculated the Gibbs free energy of coupled reactions under a wide range of substrate concentration conditions: reaction 1–5 (Table S1, Fig. S4). The best conditions suggest that the degradation of benzoate (the central metabolic intermediate in aromatic hydrocarbon degradation pathways) potentially occurs under the following metabolic conditions (a, b, and c): (a) benzoate mineralization to CO2 and H2O coupled with nitrite reduction to ammonia (ΔrG′ = −1206.3 [kJ/mol]); (b) benzoate mineralization to CO2 and H2O coupled with sulfite reduction to hydrogen sulfide (ΔrG′ = −373.6 [kJ/mol]); and (c) benzoate mineralization to CO2 and H2O coupled with nitrate reduction to nitrite (ΔrG′ = −119.9 [kJ/mol]). Coupling with sulfate reduction is not thermodynamically favorable (Table S1, Figs. S4 and S5).

The type and complexity of the exchanged substrates are another key factor that may shape the syntrophic relationship and the identity of the syntrophic partners of Lokiarchaeota. The presence of genes encoding for membrane bound electron bifurcating classes of [NiFe] hydrogenases, groups 3b and 3d, which couple the oxidation of NADH+ and NADPH+ with H2 evolution (Fig. S6) [56] also suggest there is syntrophic exchange of hydrogen between Asgard and their partners. In addition, we located genes encoding for β-oxidation enzymes (enoyl-CoA hydratase, acyl CoA dehydrogenase, and acetyl CoA acetyl transferase) and various fermentation pathways (acetate and formate) in Lokiarchaeota CR_06, which suggests that short chain fatty acid and different fermentation products could also be syntrophically exchanged (Table S4). The diverse nature of substrates could facilitate the interactions between a broader range of partners of diverse metabolic capabilities and support the conclusions driven from our thermodynamic calculations.

It is worth noting that efficient substrate and electron exchange between syntrophic partners require the presence of either biological conduits (e.g., type IV pili or flagella) or some sort of electron shuttle allowing extracellular electron transfers (e.g., multiheme cytochromes) [57, 58]. Hence, we screened CR_06 for these mechanisms and identified two candidate mechanisms for interspecies substrate and electron exchange. First, flagellar proteins are encoded by CR_06, suggesting flagella may be a potential structure mediating interspecies interactions. Second, CR_06 harbors a gene encoding for an oxidoreductase belonging to electron transfer flavoprotein-quinone oxidoreductase (CR_06_contig-100_4953_2), ETF-QO/FixC family, which potentially mediates the transfer of electrons across membranes.

Other metabolic features of Lokiarchaeota (CR_06)

The genomic analysis of Lokiarchaeota (CR_06) MAG also suggests versatile catabolic capacities potentially targeting detrital proteins and short chain fatty acids (e.g., propanoate, butyrate), which are abundant in benthic marine sediments [59, 60]. CR_06 MAG has a relatively large number of peptidases encoding genes (92 peptidases/1 Mbp) with diverse catalytic residues (e.g., aspartic, metallo, serine, etc.), which potentially degrade detrital proteins (Fig. S7). It also contains genes encoding for various classes of facilitated and active transporters, which are dedicated to shuttle oligo/dipeptides and single amino acids (e.g., polar and branched chain amino acids) across the cell membrane. Also, CR_06 encoded for enzymes enabling the utilization of wide range of amino acids (e.g., aspartate, threonine, alanine, glycine, serine, cysteine, histidine, glutamine, and lysine), which are then channeled to the central metabolic pathways and ultimately produce energy via fermentation (Fig. 3a).

MAG CR_06 suggested the capacity to breakdown short chain fatty acids e.g., propanoate, oxobutanoate as other potential substrates. Both propanoate and oxobutanoate are converted to propanoyl CoA via formate acetyl transferase and acetyl synthase, respectively. Then, the resulting propanoyl CoA is converted to acetyl-CoA via the malonyl-CoA pathway. CR_06 as well as other Asgard phyla showed different autotrophic capacities enabling carbon fixation using both the Wood–Ljungdahl (WL) pathway and the reverse tricaboxylic acid (rTCA) cycle (Fig. 3a).

Metabolic features of CR_Bathyarchaeota

The other main lineage of Archaea in these sediments is Bathyarchaeota. Bathyarchaeota CR_14 is described in detail since it has the highest quality at >95% complete and 4% contaminated. Phylogenomic analysis showed that CR_14 is clustered together with other Bathyarcaheota MAGs from CR in two distinct clades within the Bathyarchaeota phylum (Fig. 1). Also, relatively low similarity scores were observed between CR_14 and other reference Bathyarchaeota MAGs, which suggest that CR_14 belongs to a novel Bathyarchaeota class (Table S2).

Metabolic analysis revealed that CR_14 harbors genes encoding for incomplete methylotrophic methanogenesis pathway, yet is missing the key genes encoding for MCR complex (mcrABCDG). The absence of this complex suggests that CR_14 is incapable of methanogenesis and methane metabolic genes may be present to perform different functions, where they recycle methyl groups from different methylated compounds and replace the functions of methyl branch in Wood–Ljungdahl pathway (Fig. 4a). CR_14 showed the potential capability to use formate as electron and hydrogen donor through using Group 4 hydrogenases (formate hydrogenlyases) (Fig. S6) [56], and these electrons reduce trimethylamine N-oxide (TMAO) to trimethylamine via trimethylamine-oxide reductase or anaerobic dimethyl sulfoxide reductase (TMAO/DMSO reductase) [61]. Together the presence of genes encoding for trimethylamine-specific corrinoid protein as well as diverse classes of methyltransferases, CoB—CoM heterodisulfide reductase/F420 nonreducing hydrogenase (hdrABCD and mvhADG) suggest the capability of CR_14 to recycle coenzyme M (CoM) and coenzyme B (CoB) and transfer the methyl group from trimethylamine to CoM-SH [19, 62]. We located genes encoding the tetrahydromethanopterin S-methyltransferase (mtrA-H), suggesting that mtrA protein transfers the methyl group from CoM to 5,6,7,8-tetrahydromethanopterin (H4MPT) and assimilates the methyl group into acetyl-CoA via the beta subunit of CODH/ACS complex, replacing the function of the methyl branch of Wood–Ljungdahl pathway. This agrees with our finding that only the genes encoding for the carbonyl branch of Wood–Ljungdahl coupled with acetate fermentation genes (acetogenesis) were present in CR_14 and the genes encoding for methyl branch were completely missing for the same pathway. Comparative genomic analyses suggested the presence of other Bathyarchaeota lineages that are capable of coupling methylotrophy with acetogenesis (e.g., B26_1 from marine subsurface sediments and RBG 13_38_9 from Rifle sediments) (Fig. S8). Yet, CR_14 was the only Bathyarchaeota MAG having the capability to transfer methyl groups from mono-, di-, trimethylated compounds. These collective metabolic features in Bathyarchaeota CR_14 suggest this genome may represent a lineage that bridges the gap between methanogenic and acetogenic Bathyarchaeota through adopting a noncanonical acetogenic lifestyle (Table S5) [19, 20].

Fig. 4: Bathyarchaeotal metabolism and acetyl CoA synthase phylogeny.
figure 4

a Metabolic reconstruction of the Bathyarchaeota bin CR_14. Central metabolic pathways found in the genome (glycolysis, carbonyl branch of Wood–Ljungdahl, and methanogenesis-related genes) are shown in gray boxes, carbon fixation pathways (ACS/WL, PPP, and rTCA cycles) are shown in pink, electron transport chain (ETC) proteins are shown in yellow, fermentation products are shown in blue boxes, amino acids are shown in brown boxes, enzymes and enzyme complexes are shown in orange circles, energy carriers are shown in red, and metabolite and amino acid transporters are shown in light green. b Maximum-likelihood tree of the acetyl CoA synthase β subunit (ACSβ/CdhC). The tree was calculated using the best fit substitution model (LG + R9) that describes the evolutionary relationships between ACS families. The tree was made using reference sequences under the KEGG entry (K00193) collected from AnnoTree [44] and branch location was tested using 1000 ultrafast bootstraps and approximate Bayesian computation, branches with bootstrap support >80% were marked by black circles. Blue and red clades highlight sequences belong to bacterial (CdhC) and archaeal (ACSβ) versions, respectively. Scale bar indicates substitutions per site. Sequences from CR_Lokiarchaeota and CR_Bathyarchaeota bins were shaded with yellow. The tree is available with full bootstrap values in Newick format in the Supplementary Data.

Although the proposed scenario where CR_14 performs acetogenesis instead of methanogenesis is metabolically feasible, it is not clear why CR_14 invests a large amount of energy to maintain genes for methane metabolism. A plausible explanation for the presence of methanogenesis-related genes is the lack of dedicated methyltransferases that could transport methyl groups directly between methylated compounds (e.g., trimethylamine) to tetrahydromethanopterin (H4MPT). Therefore, encoding genes mediating the synthesis and cycling of CoM is necessary to use it as intermediate carrier to transport methyl groups from/to H4MPT. Also, the phylogenetic tree of the acetyl-CoA synthase beta subunit (AcsB) showed that CR_14 AcsB genes clustered together with genes recovered from other Bathyarchaeota lineages, Asgard phyla, Chloroflexi, and Altiarchaeales (Fig. 4b). The uniqueness of that clade stemmed from the previous hypothesis that Altiarchaeales members possess an ancient version of AcsB [63]. Likewise, members of Chloroflexi encoded for this archaeal version of ACSβ subunit similar to the ones present in Bathyarchaeota, Asgard phyla, and Altiarchaeales [64, 65]. Finally, this ACS type supports the hydrogen dependent, autotrophic lifestyle of the Asgard archaea [15]. Accordingly, we hypothesize that Bathyarchaeota possess the type of ACS capable of incorporating methyl groups into acetyl-CoA utilizing different carrier proteins (CoM, H4MPT, and may be other unknown carriers), allowing Bathyarchaeota to assimilate methyl groups originated from various methylated compounds.

Discussion

In this study, we employed a metagenomics-enabled genomics approach to recover MAGs from abundant archaeal phyla inhabiting deep sediments of CR margin and elucidate their potential ecological roles. A total of 31 MAGs belonging to archaeal phyla (Lokiarchaeota, Thorarchaeota, Heimdallarchaeota, Bathyarcheota, Thermoplasmatales, and Hadesarchaea) were successfully recovered from the metagenomic datasets of five different samples. Only 11 MAGs met our completion and contamination thresholds (>60% complete and <10% contamination) and were included for the detailed genomic analyses and metabolic reconstruction. More than 90% of the high-quality genomes were affiliated to Bathyarchaeota and Asgard phyla and the phylogenetic affiliations of the predicted proteins in each MAG confirmed the proper quality of the MAGs considered in this study. Remarkably, all the CR archaeal MAGs were enriched with ESP encoding genes. The wide distribution of these eukaryotic homologs indicates that ESPs are more ubiquitous in anaerobic archaea than previously recognized.

Notably, the sediments used in this study were collected from much deeper depths compared with the sediments where Bathyarchaeota and Lokiarchaeota genomes were previously reported. As such, our analyses were focused on the MAGs belonging to Bathyarchaeota and Lokiarchaeota to understand their ecological potentials under these deep marine sediment conditions.

Community interaction and Lokiarchaeota metabolic interdependencies

Deep marine sediments are rich in aliphatic and aromatic hydrocarbons that provide the associated microbes with a significant portion of their energy and carbon needs [59, 66]. Previous studies showed that hydrocarbon degradation is restricted to limited bacterial and archaeal phyla (e.g., Aminicenantes, TA06, Aerophobetes, Atribacteria, Helarchaeota, and Bathyarchaeota) [17, 19, 67]. However, the full spectrum of archaea involved in the hydrocarbon degradation processes and the nature of their interactions are not fully elucidated. In the light of the above evidence, this study provides new views regarding the ecological roles and potential metabolic capacities of Lokiarchaeota. Lokiarchaeota CR_6 expanded the range of metabolic features encoded by this phylum and predicts metabolic functions enabling the utilization of aliphatic and aromatic hydrocarbons as carbon and energy sources. Considering the high energy demands for aromatic hydrocarbon breakdown under the energy limited conditions of subseafloor sediments [54], it was unexpected to find the ATP-dependent BCR complex (class Bzd) in CR_Asgard/CR_Lokiarchaeota genomes while they have fermentation and/or acetogenic lifestyles. This strongly suggests that acetogenesis or acetate fermentation cannot be the ultimate fate of the aromatic hydrocarbon degradation process due to the low energy yields of these pathways. Under these circumstances, we propose that the CR_Lokiarchaeota members have the capacity of completely mineralizing aromatic hydrocarbons in syntrophy with microbes capable of using nitrite, sulfite, and nitrate as electron sink under the subsurface settings to increase their energy budgets and sustain their energy requirement. After testing all possible partnership scenarios based on the presence/absence profiles of the candidate pathways and the thermodynamic feasibility (Fig. 5a, b), we found the potential partners for Lokiarchaeota belonging to diverse metabolic groups (e.g., nitrate reducers, nitrite reducers, sulfate reducers, sulfite reducers, and thiosulfate reducers), however only sulfite, nitrate, and nitrite reducers are thermodynamically favored under the subsurface conditions. In the same context, a diverse group of syntrophs might require sharing substrates of different qualities (e.g., acetate, propanoate, and other short chain fatty acids). Accordingly, this suggests the presence of Lokiarchaeota syntrophic partners capable of acetate oxidation and short chain fatty acid oxidation in addition to the hydrogenotrophic ones as previously proposed [15]. Therefore, we hypothesize that thermodynamic favorability and potential diversity of the shared metabolites will push the Lokiarchaeota to syntrophy, beyond individual metabolite exchange.

Fig. 5: Overall metabolism within the environment.
figure 5

a Heat map depicts the distribution patterns of the genes encoding for the key metabolic pathways across all sampled metagenomic datasets. Metagenomes were clustered based on presence (black)/absence (white) profiles using Euclidean distance and average linkage method. Y-axis represents pathways included in the analysis and X-axis includes metagenomic datasets analyzed. b Hypothetical models for the microbial interactions in the Costa Rica sediments. The potential substrates degraded by Lokiarchaeota, Bathyarchaeota, and other microbes in the environment were colored in gray, potential metabolic products (in blue), potential interactions between Lokiarchaeota and their syntrophic partners and possible shared metabolites.

Bathyarchaeota CR_14 bridges the gap between acetogenic and methanogenic lineages

Bathyarchaeota have been shown to have both methanogenic and acetogenic capabilities [19, 20], yet the evolutionary trajectory and the ecological context driving the switch between lifestyles is not fully understood. Generally, Bathyarchaeota lineages relative to those of conventional methanogens are widely distributed in different marine benthic habitats [18], which could be facilitated by their acetogenesis capacity that degrades a wide range of organic compounds under thermodynamic favorable conditions [47].

Interestingly, phylogenomic analysis revealed that Bathyarchaeota MAGs from CR formed a novel class that are significantly divergent from previously reported acetogenic and methanogenic lineages. In these CR sediments, distinct environmental conditions may have captured the intermediate transition between acetogenic and methanogenic lifestyles. We propose that these environmental conditions include the presence of high levels of methylated compounds (e.g., methylated amines) in deep-sea environments [68] and the absence of dedicated pathways/carriers necessary to recycle these methylated compounds and directly shuttle the extracted methyl groups to the Wood–Ljungdahl pathway. The redox potential and availability of oxidized substrates in CR sediments may not favor the occurrence of full methanogenesis and the associated methanogenic archaea, as indicated by the absence of mcrA genes at the sampled depths (Fig. 5a).

In summary, this study presents archaeal MAGs recovered from deep subsurface sediments, some of which are of high quality and completeness. Numerous CR archaeal MAGs particularly those belonging to Asgard host eukaryotic signatures, yet, only the Asgard MAGs have homologs for cell division and cytoskeleton. We also report a metabolically novel Lokiarchaeota lineage capable of aliphatic and aromatic hydrocarbon degradation with a putative partnership with metabolically diverse syntrophic organisms. Also, we revealed the presence of an intermediate stage between acetogenic and methanogenic Bathyarchaeota that could convert methylated amines to acetate through linking methylotrophy to acetogenesis.