Asgard archaea capable of anaerobic hydrocarbon cycling

1Short-chain alkanes, such as methane and butane, are abundant in marine sediments and play an important role in carbon cycling with methane concentrations of ~1 Gt being processed globally through anoxic microbial communities–. Until recently, archaeal methane cycling was thought to be limited to Euryarchaeota. However, additional archaeal phyla, including Bathyarchaeota and Verstraetarchaeota, have been shown to contain proteins with homology to the activating enzyme methyl-coenzyme M reductase (Mcr) and corresponding pathways for methane utilization. Furthermore, lineages within the Euryarchaeota belonging to Candidatus Syntrophoarchaeum spp., have been shown to use methyl-CoM reductase-like enzymes for anaerobic butane oxidation. Similar to methane oxidation in many ANME-1 archaea, butane oxidation in Syntrophoarchaeum is proposed to be enabled through a syntrophic interaction with sulfur-reducing bacteria. Metagenomic reconstructions of genomes recovered from deep-sea sediments from near 2000 m depth in Guaymas Basin (GB) in the Gulf of California have revealed the presence of additional uncharacterized alkyl methyl-CoM reductase-like enzymes in metagenome-assembled genomes within the Methanosarcinales (Gom-Arc1). GB is characterized by hydrothermal alterations that transform large amounts of organic carbon into methane, polycyclic aromatic hydrocarbons, low-molecular weight alkanes and organic acids allowing for diverse microbial communities to thrive (Supplementary Table 1)–.

2Recently, genomes of a clade of uncultured archaea, referred to as the Asgard superphylum that includes the closest archaeal relatives of eukaryotes, have been recovered from anoxic environments around the world–. Diversity surveys in anoxic marine sediments show that Asgard archaea appear to be globally distributed,–. Based on phylogenomic analyses, Asgard archaea have been divided into four distinct phyla: Lokiarchaeota, Thorarchaeota, Odinarchaeota, and Heimdallarchaeota, with the latter possibly representing the closest relatives of eukaryotes. Supporting their close relationship to eukaryotes, Asgard archaea possess a wide repertoire of proteins previously thought to be unique to eukaryotes known as eukaryotic signature proteins (ESPs). These ESPs include homologs of eukaryotic proteins, which in eukaryotes are involved in ubiquitin-mediated protein recycling, vesicle formation and trafficking, endosomal sorting complexes required for transport-mediated multivesicular body formation, as well as cytokinetic abscission and cytoskeleton formation. Asgard archaea have been suggested to possess heterotrophic lifestyles and are proposed to play a role in carbon degradation in sediments; however, several members of the Asgard archaea also have genes that code for a complete Wood–Ljungdahl pathway and are therefore interesting with regard to carbon cycling in sediments,.

3Here, we present metagenome-assembled genomes (MAGs), recovered from GB deep-sea hydrothermal sediments, which represent an undescribed Asgard phylum with the metabolic potential to perform anaerobic hydrocarbon degradation using a methyl-CoM reductase-like homolog.

4We recently obtained more than ~280 gigabases of sequencing data from 11 samples taken from various sites and depths at GB hydrothermal vent sediments. De novo assembly and binning of metagenomic contigs resulted in the reconstruction of over 550 genomes (>50% complete). these genomes we detected a surprising diversity of archaea, including >20 phyla, which appear to represent up to 50% of the total microbial community in some of these samples. A preliminary phylogeny of the dataset using 37 concatenated ribosomal proteins revealed two draft genomic bins representing a previously unknown lineage of the Asgard archaea. These draft genomes, referred to as Hel_GB_A and Hel_GB_B, were re-assembled and re-binned resulting in final bins that were 82% and 87% complete and had a bin size of 3.54 and 3.84 Mbp, respectively (Table 1). An in-depth phylogenetic analysis consisting of 56 concatenated ribosomal proteins was used to confirm the placement of these final bins form a distant sister group with the Lokiarchaeota (Fig. 1a). Hel_GB_A percent abundance ranged from 3.41 × 10−3% to 8.59 × 10−5%, and relative abundance from 8.43 to 0.212. Hel_GB_B percent abundance ranged from 1.20 × 10−3% to 7.99 × 10−5%, and relative abundance from 3.41 to 0.22 compared to the total raw reads. For both Hel_GB_A and Hel_GB_B the highest abundance was seen at the site from which the bins were recovered. These numbers are comparable to other Asgard archaea whose genomes have been isolated form these sites. Hel_GB_A and Hel_GB_B had a mean GC content of 35.4% and 28%, respectively, and were recovered from two distinct environmental samples, which share similar methane-supersaturated and strongly reducing geochemical conditions (concentrations of methane ranging from 2.3 to 3 mM, dissolved inorganic carbon ranging from 10.2 to 16.6 mM, sulfate near 21 mM and sulfide near 2 mM; Supplementary Table 1) but differed in temperature (28 and 10 °C, respectively, Supplementary Table 1).

5Phylogenetic analyses of a 16S rRNA gene sequence (1058 bp in length) belonging to Hel_GB_A confirmed that it is related to Lokiarchaeota and Thorarchaeota, but is phylogenetically distinct from either of these lineages (Fig. 1b). A comparison to published Asgard archaeal 16S rRNA gene sequences indicate a phylum level division between the Hel_GB_A sequence and other Asgard archaea with a percent identity of 82.67% when compared to Lokiarchaeum GC14_75 (Supplementary Table 2). A search for ESPs in both bins revealed that they contained a similar suite of homologs compared to those previously identified in Lokiarchaeota, which is consistent with their phylogenomic relationship (Fig. 2). Yet these lineages are relatively distantly related as evidenced by their difference in GC content and relatively low-pairwise sequence identity of proteins. An analysis of the average amino acid identity (AAI) showed that Hel_GB_A and Hel_GB_B shared 1477 genes and AAI of 51.96%. When compared to Lokiarcheota_CR4, Hel_GB_A shares 634 orthologous genes out of 3595 and Hel_GB_B shares 624 out of 3157. Helarchaeota bins showed the highest AAI similarity to Odinarchaeota LCB_4 (45.9%); however, it contained fewer orthologous genes (574 out of 3595 and 555 out of 3157 for Hel_GB_A and Hel_GB_B, respectively). Additionally, the Hel_GB bins differed from Lokiarchaeota in their total gene number, for example Hel_GB_A possessed 3595 genes and CR_4 possessed 4218; this difference is consistent with the larger estimated genome size for Lokiarchaeum CR_4 compared to Hel_GB_A (~5.2 to ~4.6 Mbp) (Supplementary Table 3). These results add support to the phylum level distinction observed for Hel_GB_A and Hel_GB_B in both the ribosomal protein and 16S rRNA phylogenetic trees. We propose the name Helarchaeota after Hel, the Norse goddess of the underworld and Loki’s daughter for this lineage.

6To reconstruct the metabolic potential of these archaea, the Helarchaeota proteomes were compared to several functional protein databases (Fig. 3a). Like many archaea in marine sediments, Helarchaeota may be able to utilize organic carbon as they possess a variety of extracellular peptidases and carbohydrate degradation enzymes that include the β-glucosidase, α-l-arabinofuranosidase and putative rhamnosidase, among others (Supplementary Table 4 and Supplementary Data 1). Degraded organic substrates can then be metabolized via glycolysis and an incomplete TCA cycle from citrate to malate and a partial gamma-aminobutyric acid shunt (Fig. 3a, Supplementary Data 1). Both Helarchaeota bins are missing fructose-1,6-bisphosphatase and have few genes coding for the pentose phosphate pathway. Genes encoding for the bifunctional enzyme 3-hexulose-6-phosphate synthase/6-phospho-3-hexuloisomerase (hps-phi) were identified in Hel_GB_B suggesting they may be using the ribulose monophosphate pathway for formaldehyde anabolism. Genes coding for acetate-CoA ligase (both APM and ADP-forming) and an alcohol dehydrogenase (adhE) were identified in both genomes suggesting that the organisms may be capable of both fermentation and production of acetyl-CoA using acetate and alcohols (Supplementary Data 1). Like in Thorarchaeota and Lokiarchaeota, these genomes possess the large subunit of type IV ribulose bisphosphate carboxylase,. In addition, the Helarchaeota genomes encode for the catalytic subunit of the methanogenic type III ribulose bisphosphate carboxylase used for C-fixation. Helarchaeota are metabolically distinct from Lokiarchaeota as both Hel_GB draft genomes appear to lack a complete TCA cycle as genes coding for citrate synthase and malate/lactate dehydrogenase are absent. Both genomes also likely produce acetyl-CoA using glyceraldehyde 3-phosphate dehydrogenase which is absent in Lokiarchaeota (Supplementary Data 1). Helarchaeota genomes lack genes that code for enzymes involved in dissimilatory nitrogen and sulfur metabolism. Assimilatory genes, including sat, cysN, and cysC were found in Hel_GB_B, however, these genes were not identified in Hel_GB_A. This absence may be indicative of species-specific characteristics or could be a results of genome incompleteness. Additional genomes of members of the Helarchaeota will help to fully understand the diversity of these pathways across the whole phylum.

7Interestingly, both Helarchaeota genomes have mcrABG-containing gene clusters encoding putative methyl-CoM reductase-like enzymes (Fig. 3b, Supplementary Figure 1),,. Phylogenetic analyses of both the A subunit of methyl-CoM reductase-like enzymes (Supplementary Figure 2) as well as the concatenated A and B subunits (Fig. 3b) revealed that the Helarchaeota sequences are distinct from those involved in methanogenesis and methane oxidation but cluster with homologs from butane-oxidizing Syntrophoarchaea and Bathyarchaeota with high-statistical support (rapid bootstrap support/single-branch test bootstrap support/posterior probability of 99.8/100/1; Fig. 3b) excluding the distant homolog of Ca. Syntrophoarchaeum caldarius (OFV68676). Analysis of the Helarchaeota mcrA alignment confirmed that amino acids present at their active sites are similar to those identified in Bathyarchaeota and Syntrophoarchaeum methyl-CoM reductase-like enzymes (Supplementary Figure 3). In Syntrophoarchaeum, the methyl-CoM reductase-like enzymes have been suggested to activate butane to butyl-CoM. It is proposed that this process is then followed by the conversion of butyl-CoM to butyryl-CoA; however, the mechanism of this reaction is still unknown. Butyryl-CoA can then be oxidized to acetyl-CoA that can be further feed into the Wood–Ljungdahl pathway to produce CO2. While some n-butane is detected in GB sediments (usually below 10 µM), methane is the most abundant hydrocarbon (Supplementary Table 1) followed by ethane and propane (often reaching the 100 µM range); thus, a spectrum of short-chain alkanes could potentially be metabolized by Helarchaeota.

8Next, we searched for genes encoding enzymes potentially involved in hydrocarbon utilization pathways, including propane and butane oxidation. Along with the methyl-CoM reductase-like enzyme that could convert alkane to alkyl-CoM, Helarchaeota possess heterodisulfide reductase subunits ABC (hdrABC), which is needed to recycle the CoM and CoB heterodisulfides after this reaction occurs (Figs. 3 and 4),. The conversion of alkyl-CoM to acyl-CoA is currently not understood in archaea capable of butane oxidation. Specialized alkyl-binding versions of methyltransferases would be required to convert alkyl-CoM to butyl-CoA or other acyl-CoAs, as discussed for Ca. S. butanivorans. Genes coding for methyltransferases were identified in both Helarchaeota genomes, including a likely tetrahydromethanopterin S-methyltransferase subunit H (MtrH) homolog (Fig. 4; Supplementary Data 1). Short-chain acyl-CoA could be oxidized to acetyl-CoA using the beta-oxidation pathway via a short-chain acyl-CoA dehydrogenase, enoyl-CoA hydratase, 3-hydroxyacyl-CoA dehydrogenase, and acetyl-CoA acetyltransferase, candidate enzymes for all of which are present in the Helarchaeota genomes and are also found in genomes of other Asgard archaea (Fig. 4). Along with these enzymes, genes coding for the associated electron transfer systems, including an Fe–S oxidoreductase and all subunits of the electron transfer flavoprotein complex were identified in Helarchaeota (Fig. 4). Acetyl-CoA produced by beta-oxidation might be further oxidized to CO2 via the Wood–Ljungdahl pathway, using among others the classical 5,10-methylene-tetrahydromethanopterin reductase (Figs. 3a and 4).

9To make anaerobic alkane oxidation energetically favorable, it must be coupled to the reduction of an internal electron acceptor or transferred to a syntrophic partner that can perform this reaction,,. We could not identify an internal electron sink or any canonical terminal reductases used by ANME archaea (such as iron, sulfur, or nitrogen reductases), leading to the conclusion that a syntrophic partner organism would be necessary to enable growth on short-chain hydrocarbons. However, we could not identify any obvious syntrophic partner organisms based on co-occurrence analyses of abundance profiles in metagenomic datasets generated in this study.

10An evaluation of traditional energy transferring mechanisms showed that the Helarchaeota bins lack genes coding for NADH:ubiquinone oxidoreductase, F420-dependent oxidoreductase, F420H2:quinone oxidoreductase and NADH:quinone oxidoreductase that were identified in Ca. S. butanivorans (Fig. 4). These protein complexes are important for energy transfer across the cell membrane and are common among syntrophic organisms,,. Helarchaeota also lack genes coding for pili or cytochromes that are often involved in direct electron transfer to a bacterial partner, as demonstrated for different ANME archaea,. Therefore, Helarchaeota may use a thus far unknown mechanism for energy conservation. Below we analyzed potential energy-transferring mechanisms that might be involved in syntrophic interactions between Helarchaeota and potential partner organisms.

11A possible candidate for energy transfer to a partner may be formate dehydrogenase because substrate exchange in form of formate has previously been described to occur between methanogens and sulfur-reducing bacteria. Helarchaeota genomes code for the alpha and beta subunits of a membrane-bound formate dehydrogenase (EC. 1.2.1.2) that could facilitate this transfer (Fig. 2, Supplementary Data 1). However, to our knowledge formate transfer has not been shown to mediate methane oxidation. Alternatively, Helarchaeota may possess a previously undiscovered redox-active complex. In both Helarchaeota bins, a gene cluster was found encoding three proteins that were identified as members of the HydB/Nqo4-like superfamily, Oxidored_q6 superfamily, and a Fe–S disulfide reductase with a FlpD domain (mvhD) (Fig. 5a). An analysis of these three proteins showed that each possessed transmembrane motifs (Fig. 5b, and Supplementary Methods). While the membrane association of the disulfide reductase/FlpD needs to be confirmed, interactions with the other two membrane-associated subunits may allow for the bifurcated electrons to be transferred across the membrane.

12Finally, hydrogen production and release was also considered as a possible electron sink for Helarchaeota. We identified several hydrogenase subunits and putative Fe-S disulfide reductase-encoding genes in the Helarchaeota genomes. Subsequent phylogenetic analyses revealed that the majority of these hydrogenases represent small and large subunits of group IIIC hydrogenases (methanogenic F420-non-reducing hydrogenase (mvh)) that are usually involved in bifurcating electrons from hydrogen (Supplementary Figure 4, Supplementary Data 1). In contrast, while homologs belonging to the above mentioned Oxidored_q6 superfamily protein family are often found to be associated with group IV hydrogenases, canonical membrane-bound group IV-hydrogenases could not be identified in the genomes of the Helarchaeota. Altogether, this indicates that hydrogen could play a central role in energy metabolism of Helarcharota, but the absence of a classical membrane-bound hydrogenase makes it unlikely that hydrogen is the major syntrophic electron carrier.

13Historically methanogenesis and anaerobic methane oxidation were regarded as the only examples of anaerobic archaeal short-chain alkane metabolism. The enzymes acting in these pathways were considered to be biochemically and phylogenetically unique and limited to lineages within the Euryarchaeota. This study represents the discovery of the previously unknown phylum referred to as Helarchaeota, whose members encode a mcr-like gene cluster. This opens the possibility that some representatives of the Asgard archaea may have the potential for anaerobic short-chain alkane oxidation. Since the presence of these mcr genes is restricted to Helarchaeota among the known Asgard archaea, these genes were likely transferred to Helarchaeota and do not constitute an ancestral trait within the Asgard superphylum. Based on current phylogenetic analysis, the Helarchaeota mcr gene cluster may have been horizontally acquired from either Bathyarchaeota or Ca. Syntrophoarchaeum (Fig. 1b, Supplementary Figure 3). Due to this close relationship, we based our analysis of Helarchaeota’s ability to perform anaerobic short-chain hydrocarbon oxidation on the pathway proposed for Ca. Syntrophoarchaeum. Helarchaeota probably utilize a similar short-chain alkane as a substrate in lieu of methane, but given the low-butane concentrations at our site it may not be the only substrate.

14Our comparison to Ca. S. butanivorans shows a consistent presence in genes necessary for this metabolism including a complete Wood–Ljungdahl pathway, acyl oxidation pathway, and internal electron transferring systems. Some of these electron-transferring systems are essential housekeeping components that may act as electron carriers for oxidation reactions. Interestingly, in the Wood–Ljungdahl pathway identified in Ca. S. butanivorans, the bacterial enzyme 5,10-methylene-tetrahydrofolate reductase (met) is thought to be substituting for the missing 5,10-methylene-tetrahydromethanopterin reductase (mer). In contrast, Helarchaeota encode the canonical archaeal-type mer. To render anaerobic butane oxidation energetically favorable, it must be coupled to the reduction of an electron acceptor such as nitrate, sulfate or iron,,. In ANME archaeum that lack genes for internal electron acceptors, methane oxidation is enabled through the transfer of electrons to a syntrophic partner organism. In Syntrophoarchaeum, syntrophic butane oxidation is thought to occur through the exchange of electrons via pili and/or cytochromes with sulfate-reducing bacteria. Helarchaeota do not appear to encode any of the systems traditionally associated with syntrophy and no partner was identified in this study. Thus, further research is needed to identify possible bacterial partners.

15Furthermore, the hypothesis that Helarchaeota have the ability to utilize short-chain alkanes remains to be confirmed as the genomes of members of this group do not encode canonical routes for electron transfer to a partner bacterium. However, we identified potential enzymes that may be involved in transfer of electrons. Some methanogenic archaea use formate for syntrophic energy transfer to a syntrophic partner; therefore, the reverse reaction has been speculated to be energetically feasible for methane oxidation. If this is true, the presence of a membrane-bound formate dehydrogenase in the Helarchaeota genomes may support this electron-transferring mechanism, however, to our knowledge this has never been shown for an ANME archaea so far. Alternatively, the type 3 NiFe-hydrogenases encoded by Helarchaeota may be involved in transfer of hydrogen to a partner organism. For example, we identified a protein complex distantly related to the mvh–hdr of methanogens for electron transfer (Supplementary Methods). Mvh–hdr structures have been proposed to be potentially used by facultative hydrogenotrophic methanogens for energy transfer, but the directionality of hydrogen exchange could easily be reversed. These methanogens form syntrophic associations with fermenting, H2-producing bacteria, lack dedicated cytochromes or pili and use the mvh–hdr for electron bifurcation. The detection of a hydrophobic region in the mvh–hdr complex led to the suggestion that this complex could be membrane bound and act as mechanism for electron transfer across the membrane; however, a transmembrane association has never been successfully shown. While the membrane association of the disulfide reductase/FlpD needs to be confirmed, we were able to detect several other transmembrane motifs in the associated proteins that could potentially allow electron transfer in form of hydrogen to an external partner. Thus, while we propose that the most likely explanation for anaerobic short-chain alkane oxidation in Helarchaeota is via a syntrophic interaction with a partner, additional experiments are needed to confirm this working hypothesis.

16The discovery of alkane-oxidizing pathways and possible syntrophic interactions in a phylum of Asgard archaea indicates a much wider phylogenetic range for hydrocarbon utilization. Based on phylogenetic analyses it seems most likely that the Helarchaeota mcr operon may have been horizontally transferred from either Bathyarchaeota or Syntrophoarchaea. However, the preservation of a horizontally transferred pathway is indicative of a competitive advantage; it follows that gene transfers among different archaeal phyla reflect alkane oxidation as a desirable metabolic trait. The discovery of the alkyl-CoM reductases and alkane-oxidizing pathways among the Asgard archaea indicates ecological roles for these still cryptic organisms, and opens up a wider perspective on the evolution and expansion of hydrocarbon-oxidizing pathways throughout the archaeal domain.

17Samples analyzed here are part of a study that aimed to characterize the geochemical conditions and microbial community of GB hydrothermal vent sediments (Gulf of California, Mexico),. The two genomic bins discussed in this paper, Hel_GB_A and Hel_GB_B, were obtained from sediment core samples collected in December 2009 on Alvin dives 4569_2 and 4571_4, respectively. Immediately after the dive, freshly recovered sediment cores were separated into shallow (0–3 cm), intermediate (12–15 cm), and deep (21–24 cm) sections for further molecular and geochemical analysis, and frozen at −80 °C on the ship until shore-based DNA extraction. Hel_GB_A was recovered from the intermediate sediment (~28 °C) and Hel_GB_B was recovered from shallow sediment (~10 °C) from a nearby core (Supplementary Table 1); the sampling context and geochemical gradients of these hydrothermally influenced sediments are published and described in detail,.

18DNA was extracted from sediment samples using the MO BIO—PowerMax Soil DNA Isolation kit and sent to the Joint Genome Institute (JGI) for sequencing.

19A half a lane of Illumina reads (HiSeq-2500 1TB, read length of 2 × 151 bp) were generated at Joint Genome Institute for each sample, producing a total of 226,647,966 and 241,605,888 reads from dives 4569-2 and 4571-4, respectively. The average percent of reads with a phred-score (Q) ≥ 30 was 86.2% and 90.39% and the average base quality score was 34.35 ± 7.73 and 35.38 ± 6.52 for samples from dive 4569-2 and 4571-4, respectively. The JGI performed read quality checks and generated a first assembly using the following methods: BBDuk adapter trimming removed known Illumina adapters. The reads were further processed using BBDuk quality filtering and trimming to remove reads with a quality score less than 12, containing more than three “Ns”, or with quality scores (before trimming) averaging less than 3 over the read length, or length under 51 bp after trimming. In addition, reads matching Illumina artifacts or phiX were discarded. The remaining reads were mapped to a masked version of the human HG19 with BBMap and all hits over 93% sequence identity to the human genome were discarded. Trimmed, screened, paired-end Illumina reads were assembled using the megahit assembler using a range of kmers. Assemblies were preformed with default parameters in megahit with the following options: “–k-list 23,43,63,83,103,123”. High-quality reads were mapped to the final assembly to calculate coverage information using bbmap by excluding all parameters except ambiguous=random as described by JGI.

20The contigs from the JGI assembled data were binned using ESOM, MetaBAT, and CONCOCT and resulting bins were combined using DAS Tool (version 1.0). For ESOM, binning was performed on contigs with a minimum length of 2000 bp using the K-batch algorithm for training after running the perl script esomWrapper.pl. Emerging self-organizing maps (ESOM) were manually sorted and curated. The bins were extracted using getClassFasta.pl (using −loyal 51). Reference genomes were included to add genetic signatures for the assembled contigs and improve binning. For CONCOCT, Anvi’o (v2.2.2) was used as the metagenomic workflow pipeline. Coverage information was obtained by mapping all high-quality reads of each sample against the assembly of another sample using the BWA-MEM algorithm in paired-end mode (bwa-0.7.12-r1034; using default settings). The resulting sam file was sorted and converted to bam using samtools (version 0.1.19). The bam file was prepared for Anvi’o using the script anvi-init-bam and a contigs database generated using anvi-gen-contigs-database. These files were the input for anvi-profile. Generated profiles for the assemblies were combined using anvi-merge and the resulting bins summarized using anvi-summarize (-C CONCOT). If not mentioned otherwise, the scripts were used with default settings. Metabat was also used as a binning approach (v1). As described for Anvi’o the input consisted of the scaffold files (≥2000 bp) and the mapping files. First, each of the mapping files were summarized using jgi_summarize_bam_contig_depths and then metabat was run using the following settings: –minProb 75 –minContig 2000 –minContigByCorr 2000. Results from the three different binning tools were combined using DAS Tool (version 1.0). For each of the binning tools a scaffold-to-bin list was prepared and DAS Tool run on each of the eleven scaffold files as follows: DAS_Tool.sh -i Anvio_contig_list.tsv,Metabat_contig_list.tsv,ESOM_contig_list.tsv -l Anvio,Metabat,ESOM -c scaffolds.fasta –write_bins 1. CheckM lineage_wf (v1.0.5) was run on bins generated from DAS_Tool and 577 bins showed an completeness >50% and were characterized further. 37 Phylosift identified marker genes were used for preliminary phylogenetic identification of individual bins (Supplementary Table 5). Thereby, we identified two genomes, belonging to a previously uncharacterized phylum within the Asgard archaea, which we named Helarchaeota. To improve the quality of the two Helarchaeota genomes IDBA-UD was run on raw data using the command: “idba_ud -r Guay9_METAGENOME.fasta -o G9 –pre_correction –mink 75 –maxk 105 –seed_kmer 55 –num_threads 30”. Metaspades was run on Raw data and Metabat assembled bins using as follows: “metaspades.py –12 Guay16.11400.5.204846.CTCTCTA-CGTCTAA.filter-METAGENOME.fastq -o Metaspades –only-assembler –meta”. Binning procedures (using scaffolds longer than 2000 bp) as described above for the original bins were repeated with these redone assembles. All bins were compared to the original Helarchaeota bins using blastn for identification. Mmgenome and CheckM were used to calculate genome statistics (i.e., contig length, genome size, contamination, and completeness). The highest quality Helarchaeota bin from each sample was chosen for further analyses. For the 4572-4 dataset, the best bin was generated using the Metaspades reassembly on the trimmed data and for the 4569-2 dataset the best bin was recovered using the Metaspades reassembly on the original Hel bin contigs. The final genomes were further cleaned by GC content, paired-end connections, sequence depth and coverage using Mmgenome. CheckM was rerun on cleaned bins to estimate the Hel_GB_A to be 82% and Hel_GB_B to be 87% complete and both bins were characterized by a low degree of contamination (between 1.4 and 2.8% with no redundancy) (Table 1). Genome size was estimated to be 4.6 Mbp for Hel_GB_A and 4.1 for Hel_GB_B and was calculated using percent completeness and bin size to extrapolate the likely size of the complete genome. CompareM was used to analysis differences between Helarchaeota bins and published Asgard bins using the command python comparem aai_wf –tmp_dir tmp/ –file_ext fa -c 8 aai_compair_loki aai_compair_loki_output (https://github.com/dparks1134/CompareM). Read abundance summarized by jgi_summarize_bam_contig_depths were used to calculate relative read abundance and total percent of metagenomic reads. Relative read abundance was calculated as total read abundance normalized to genome size and divided by total reads. Relative read abundance was then multiplied by the constant 1 × 1012 for clarity. Total percent of metagenomic reads was calculated as total read abundance divided by total reads times 100. Relative read abundance was compared to other genomics bins recovered from these sites to look for co-occurrence.

21Neither bin possessed a 16S rRNA gene sequence, and to uncover potentially unbinned 16S rRNA gene sequences from Helarchaeota, all 16S rRNA gene sequences obtained from samples 4569_2 and 4571_4 were identified using JGI-IMG annotations, regardless of whether or not the contig was successfully binned. These 16S rRNA gene sequences were compared using blastn (blastn -outfmt 6 -query Hel_possible_16s.fasta –db Hel_16s -out Hel_possible_16s_blast.txt -evalue 1E-20) to recently acquired 16S rRNA gene sequences from MAGs recovered from preliminary data from additional GB sites. A 37 Phylosift marker genes tree was used to assign taxonomy to these MAGs. We were able to identify five MAGs that possessed 16S and that formed a monophyletic group with our Hel_GB bins (Supplementary Table 2; Megxx in Fig. 2). Of the unbinned 16S rRNA gene sequences one was identified as likely Helarchaeota sequence. The contig was retrieved from the 4572_4 assembly (designated Ga0180301_10078946) and was 2090 bp long and encoded for an 16S rRNA gene sequence that was 1058 bp long. Given the small size of this contig relative to the length of the 16S rRNA gene none of the other genes on the contig could be annotated. Blastn comparison to published Asgard 16S rRNA gene sequences was performed using the following command: blastn -outfmt 6 -query Hel_possible_16s.fasta –db Asgrad_16s -out Hel_possible_16s_blast.txt -evalue 1E-20 (Supplementary Table 2). The GC content of each 16S rRNA gene sequence was calculated using the Geo-omics script length+GC.pl (https://github.com/Geo-omics/scripts/blob/master/AssemblyTools/length%2BGC.pl). For a further phylogenetic placement, the 16S rRNA gene sequences were aligned to the SILVA database (SINA v1.2.11) using the SILVA online server and Geneious (v10.1.3) was used to manually trim sequences. The alignment also contained 16S rRNA gene sequences from the preliminary Helarchaeota bins. The cleaned alignment was used to generated a maximum-likelihood tree with RAxML as follows: “/raxmlHPC-PTHREADS-AVX -T 20 -f a -m GTRGAMMA -N autoMRE -p 12345 -x 12345 -s Nucleotide_alignment.phy -n output” (Fig. 1b).

22For a more detailed phylogenetic placement, we used BLASTp to identify orthologs of 56 ribosomal proteins in the two Helarchaeota bins, as well as from a selection of 130 representative taxa of archaeal diversity and 14 eukaryotes. The full list of marker genes selected for phylogenomic analyses is shown in Supplementary Table 6. Individual protein datasets were aligned using mafft-linsi and ambiguously aligned positions were trimmed using BMGE (-m BLOSUM30). Maximum likelihood (ML) individual phylogenies were reconstructed using IQtree v. 1.5.5 under the LG+C20+G substitution model with 1000 ultrafast bootstraps that were manually inspected. Trimmed alignments were concatenated into a supermatrix, and two additional datasets were generated by removing eukaryotic and/or DPANN homologs to test the impact of taxon sampling on phylogenetic reconstruction. For each of these concatenated datasets, phylogenies were inferred using ML and Bayesian approaches. ML phylogenies were reconstructed using IQtree under the LG+C60+F+G+PMSF model. Statistical support for branches was calculated using 100 bootstraps replicated under the same model. To test robustness of the phylogenies, the dataset was subjected to several treatments. For the “full dataset” (i.e., with all 146 taxa), we tested the impact of removing the 25% fastest-evolving sites, as within a deep phylogenetic analysis, these sites are often saturated with multiple substitutions and, as a result of model-misspecification can manifest in an artifactual signal–. The corresponding ML tree was inferred as described above. Bayesian phylogenies were reconstructed with Phylobayes for the dataset “without DPANN” under the LG+GTR model. Four independent Markov chain Monte Carlo were run for ~38,000 generations. After a burn-in of 20%, convergence was achieved for three of the chains (maxdiff < 0.29). The initial supermatrix was also recoded into four categories, in order to ameliorate effects of model misspecification and saturation and the corresponding phylogeny was reconstructed with Phylobayes, under the CAT+GTR model. Four independent Markov chain Monte Carlo chains were run for ~49,000 generations. After a burn-in of 20 convergence was achieved for all four the chains (maxdiff < 0.19). All phylogenetic analyses performed are summarized in Supplementary Table 7, including maxdiff values and statistical support for the placement of Helarchaeota, and of eukaryotes.

23McrA homologs were aligned using mafft-linsi, trimmed with trimAL, and the final alignment consisting of 528 sites was subjected to phylogenetic analyses using IQtree v. 1.5.5 with the LG+C60+R+F model. Support values were estimated using 1000 ultrafast boostraps and SH-like approximate likelihood ratio test, respectively. Sequences for McrA and B were aligned separately with mafft-linsi and trimmed using trimAL Subsequently, McrA and McrB encoded in the same gene cluster, were concatenated yielding a total alignment of 972 sites. Bayesian and ML phylogenies were inferred using IQtree v. 1.5.5 with the mixture model LG+C60+R+F and PhyloBayes v. 3.2 using the CAT-GTR model. For ML inference, support values were estimated using 1000 ultrafast boostraps and SH-like approximate likelihood ratio test, respectively. For Bayesian analyses, four chains were run in parallel, sampling every 50 points until convergence was reached (maximum difference < 0.07; mean difference < 0.002). The first 25% or the respective generations were selected as burn-in. Phylobayes posterior predictive values were mapped onto the IQtree using sumlabels from the DendroPy package. The final trees were rooted artificially between the canonical Mcr and divergent Mcr-like proteins, respectively. Original alignment and treefiles are available upon request.

24Gene prediction for the two Helarchaeota bins was performed using prodigal (V2.6.2) with default settings and Prokka (v1.12) with the extension “–kingdom archaea”. Results for both methods were comparable and yielded a total of 3574–3769 and 3164–3287 genes for Hel_GB_A and Hel_GB_B, respectively, with Prokka consistently identifying fewer genes. Genes were annotated by uploading the protein fasta files from both methods to KAAS (KEGG Automatic Annotation Server) for complete or draft genomes to assign orthologs. Files were run using the following settings: prokaryotic option, GhostX, and bi-directional best hit (BBH). Additionally, genes were annotated by JGI-IMG to confirm hits using two independent databases. Hits of interest were confirmed using blastp on the NCBI webserver. The dbCAN and MEROPS webserver were run using default conditions for identification of carbohydrate degrading enzymes and peptidases respectively. Hits with e-values lower than e−20 were discarded. In addition to these methods an extended search was used to categorize genes involved in butane metabolism, syntrophy and energy transfer.

25Identified genes predicted to code for putative alkane oxidation proteins were similar to those described from Candidatus Syntrophoarchaeum spp. Therefore, a blastp database consisting of proteins predicted to be involved in the alkane oxidation pathway of Ca. Syntrophoarchaeum was created in order to identify additional proteins in Helarchaeota, which may function in alkane oxidation. Positive hits were confirmed with blastp on the NCBI webserver and compared to the annotations from JGI-IMG, Interpro, Prokka, and KAAS annotation. Genes for mcrABG were further confirmed by a HMMER search to a published database using the designated threshold values and multiple MCR trees (see Methods). To confirm that the contigs with the mcrA gene cluster were not missbined, all other genes on these contigs were analyzed for their phylogenetic placement and gene content. The prodigal protein predictions for genes on the contigs with mcrA operons were used to determine directionality and length of the potential operon.

26To identify genes that are involved in electron and hydrogen transfer across the membrane, a database was created of known genes relevant in syntrophy that were download from NCBI. The protein sequences of the two Helarchaeota genomes were blasted against the database to detect relevant hits (E-value ≥ e−10). All hits were confirmed using the NCBI webserver, Interpro, JGI-IMG, and KEGG. Hydrogenases were identified by a HMMER search to published database using the designated threshold values. Hits were confirmed with comparisons against JGI annotations and NCBI blasts, the HydDB database and a manual database made from published sequences,. All detected hydrogenases were used to generate two phylogenetic trees, one for proteins identified as small subunits and one for large subunits in order to properly identify the different hydrogenase subgroups. Hydrogenases that are part of the proposed complex were then further analyzed to evaluate if this was a possible operon by looking for possible transcription factors and binding motifs (Supplementary Methods).

27Gene prediction for the two Helarchaeota bins was performed using prodigal (V2.6.2) with default settings. All the hypothetical proteins inferred in both Helarchaeaota were used as seeds against InterPro, arCOG, and nr using BLAST. The annotation table from Zaremba-Niedzwiedzka et al. was used as a basis for the comparison. The IPRs (or in some cases, the arCOGs) listed in the Zaremba-Niedzwiedzka et al. were searched for in the Helarchaeota genomes, and the resulting information was used to complete the presence/absence of table. When something that had previously been detected in an Asgard bin was not found in a Helarchaeota bin using the InterPro/arCOG annotations, BLASTs were carried out using the closest Asgard seeds to verify the absence. In some cases, specific analyses were used to verify the homology or relevance of particular sequences. The details for each individual ESP are depicted in Supplementary Methods.

28Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.