Metabolic Potential for Reductive Acetogenesis and a Novel Energy-Converting [NiFe] Hydrogenase in Bathyarchaeia From Termite Guts – A Genome-Centric Analysis

Symbiotic digestion of lignocellulose in the hindgut of higher termites is mediated by a diverse assemblage of bacteria and archaea. During a large-scale metagenomic study, we reconstructed 15 metagenome-assembled genomes of Bathyarchaeia that represent two distinct lineages in subgroup 6 (formerly MCG-6) unique to termite guts. One lineage (TB2; Candidatus Termitimicrobium) encodes all enzymes required for reductive acetogenesis from CO2 via an archaeal variant of the Wood–Ljungdahl pathway, involving tetrahydromethanopterin as C1 carrier and an (ADP-forming) acetyl-CoA synthase. This includes a novel 11-subunit hydrogenase, which possesses the genomic architecture of the respiratory Fpo-complex of other archaea but whose catalytic subunit is phylogenetically related to and shares the conserved [NiFe] cofactor-binding motif with [NiFe] hydrogenases of subgroup 4 g. We propose that this novel Fpo-like hydrogenase provides part of the reduced ferredoxin required for CO2 reduction and is driven by the electrochemical membrane potential generated from the ATP conserved by substrate-level phosphorylation; the other part may require the oxidation of organic electron donors, which would make members of TB2 mixotrophic acetogens. Members of the other lineage (TB1; Candidatus Termiticorpusculum) are definitely organotrophic because they consistently lack hydrogenases and/or methylene-tetrahydromethanopterin reductase, a key enzyme of the archaeal Wood–Ljungdahl pathway. Both lineages have the genomic capacity to reduce ferredoxin by oxidizing amino acids and might conduct methylotrophic acetogenesis using unidentified methylated compound(s). Our results indicate that Bathyarchaeia of subgroup 6 contribute to acetate formation in the guts of higher termites and substantiate the genomic evidence for reductive acetogenesis from organic substrates, possibly including methylated compounds, in other uncultured representatives of the phylum.


INTRODUCTION
Although Bathyarchaeia are widespread in anoxic environments, their physiology is only poorly understood. In the absence of any isolates and with only a few microscopic observations of their cells (Collins et al., 2005;Kubo et al., 2012), our knowledge about this deep-branching lineage is based almost exclusively on amplicon libraries of archaeal 16S rRNA genes and metagenomic studies (reviewed by Zhou et al., 2018).
Ribosomal RNA genes affiliated with the Miscellaneous Crenarchaeotal Group (MCG) had already been recovered in early analyses of archaeal diversity in diverse anoxic habitats (e.g., Schleper et al., 1997;Inagaki et al., 2003;Ochsenreiter et al., 2003), including the intestinal tract of termites (Friedrich et al., 2001). Meanwhile, an enormous diversity of sequences from this group, which comprises numerous deep-branching lineages, has been recovered from a wide range of marine and freshwater habitats and terrestrial environments (e.g., Kubo et al., 2012;Fillol et al., 2016). A few years ago, the MCG was elevated to the phylum level (Bathyarchaeota; Meng et al., 2014), but the most recent genome-based taxonomy demoted them again to the class level (Bathyarchaeia; Rinke et al., 2020). While the rank of the taxon is not relevant in the current context, we maintained the subgroup numbering used in previous studies (e.g., Kubo et al., 2012;Lazar et al., 2016) but replaced the prefix "MCG-" with the prefix "Bathy-" (Yu T. et al., 2018).
The abundance of Bathyarchaeia in many anoxic habitats implies potentially important roles in biogeochemical cycles (Evans et al., 2015;He et al., 2016). Reconstruction of metagenome-assembled genomes (MAGs) provided information concerning the metabolic capacities of Bathyarchaeia and inspired predictions of their putative roles in anoxic sediments (reviewed by Zhou et al., 2018). Several studies suggested that Bathyarchaeia are organotrophic and utilize a variety of organic substrates (e.g., Meng et al., 2014;He et al., 2016;Lazar et al., 2016). The discovery of genes encoding a methyl-coenzyme M reductase (Mcr) complex and a complete Wood-Ljungdahl pathway in bathyarchaeon BA1 provided the first evidence of methanogenesis outside the Euryarchaeota (Evans et al., 2015). Other studies detected key enzymes of the pathway in bathyarchaeal genomes of several subgroups and proposed that these lineages are involved in reductive acetogenesis from CO 2 (He et al., 2016;Lazar et al., 2016).
Considering the putative roles of Bathyarchaeia in methanogenesis and reductive acetogenesis and the evidence for the utilization of lignin-derived methoxy groups (Yu T. et al., 2018), the presence of this group in termite guts is intriguing. Termites efficiently digest wood and other lignocellulosic substrates, either sound or in different stages of humification (Brune, 2014), in symbiosis with a specialized gut microbiota housed in their enlarged hindgut compartments (Brune and Dietrich, 2015). Hydrogen produced in microbial fermentation processes serves as an electron donor for the reduction of CO 2 , yielding acetate and methane as major products (Breznak and Switzer, 1986;Brauman et al., 1992). Methanogenesis in termite guts involves a diverse assemblage of hydrogenotrophic and methyl-reducing archaea (Brune, 2018), but reductive acetogenesis, which can contribute up to two-thirds of total acetate production, has so far been considered a bacterial activity.
In lower termites, reductive acetogenesis has been attributed to acetogenic members of the phylum Spirochaetes) (e.g., Leadbetter et al., 1999;Ohkuma et al., 2015) and a novel lineage of uncultured Deltaproteobacteria (Rosenthal et al., 2013;Ikeda-Ohtsubo et al., 2016). In higher termites (family Termitidae), which diverged from the lower termites about 50 million years ago (Bucek et al., 2019), the situation is more complex. Particularly in the humus-feeding and soil-feeding groups, where the potential rates of reductive acetogenesis decrease in favor of methanogenesis (Brauman et al., 1992;Tholen and Brune, 1999), spirochetes are less abundant than in wood-feeding groups . A study based on the formyltetrahydrofolate synthetase (FTHFS) gene, a key enzyme of the Wood-Ljungdahl pathway that has been used as a marker for reductive acetogenesis, indicated that the community of potential acetogens shifts from spirochetes in lower termites to clostridia in higher termites (Ottesen and Leadbetter, 2011).
In a large-scale metagenomic study of the gut microbiota of eight higher termites, we obtained 15 MAGs assigned to Bathyarchaeia (Hervé et al., 2020). Preliminary analysis revealed that they fell into a cluster comprising mainly termite gut MAGs, with members of Bathy-1 and Bathy-6 as next relatives. Here, we conducted detailed phylogenomic analyses of these MAGs and investigated their potential capacity for methanogenesis and reductive acetogenesis using a genome-centric approach.

Phylogeny of Termite Gut Bathyarchaeia
Bathyarchaeal MAGs were recovered from seven of the eight higher termites investigated, regardless of their feeding group (Hervé et al., 2020; Table 1). Their absence from Microcerotermes parvus is most likely caused by the low total number of MAGs obtained from the metagenomes of this species. Based on average nucleotide identity (ANI), the MAGs were assigned to nine phylotypes ( Table 1). MAGs of the same phylotype were always derived from different gut compartments of the same host species, indicating that they most likely represent bathyarchaeal populations distributed along the entire hindgut. Eleven of the 15 MAGs fulfill the criteria for high-quality MAGs (>90% complete and <5% contamination; Bowers et al., 2017). Except for phylotype 5, each phylotype is represented by at least one high-quality MAG, which allows robust inference of metabolic potentials (Nelson et al., 2020).
Phylogenomic analysis placed all phylotypes from termite guts within subgroup Bathy-6, an apical lineage of Bathyarchaeia that is well represented mostly in 16S rRNA gene libraries (He et al., 2016) but comprises only a few MAGs from marine or estuarine sediments and the deep subsurface (Figure 1). The MAGs from termite guts form two distinct lineages, TB1 (phylotypes 1-7) and TB2 (phylotypes 8 and 9). TB2 is a sister group of bathyarchaeon SZUA-568 (hereafter denoted as Bathy-6-S), a MAG retrieved from marine hydrothermal vent sediments. Other MAGs in the radiation of Bathy-6 are a Average nucleotide identity (ANI) = 99%; for details, see Supplementary Figure S1. b The first letters of the MAG names indicate the host species (Co, Cornitermes sp.; Emb, Embiratermes neotenicus; Lab, Labiotermes labralis; Th, Termes hospes; Cu, Cubitermes ugandensis; Nc, Nasutitermes corniger; Nt, Neocapritermes taracua). c Relative abundance of the reads assigned to each MAG among the total number of reads in the corresponding metagenome (Hervé et al., 2020). d Completeness and contamination were estimated with CheckM using 107 single-copy marker genes . For detailed results of the CheckM analysis, see Supplementary Table S1. e For NCBI Nucleotide database; IMG genome IDs are given in Supplementary Table S1. f Referred to as phylotypes Bathy-6-S (J. Pan and Z. Zhou, unpublished), Bathy-6-B (Harris et al., 2018), and Bathy-6-A (Lazar et al., 2016).
FIGURE 1 | Genome-based phylogeny of termite gut Bathyarchaeia, illustrating the relationship of lineages TB1 and TB2 to other MAGs in the Bathy-6 subgroup (f__UBA233 in the GTDB taxonomy). MAGs of other subgroups that are mentioned in the text are marked in bold. The maximum-likelihood tree was inferred from a concatenated alignment of 43 marker genes using the LG+F+I+G4 model and rooted with selected Crenarchaeota and Euryarchaeota as outgroup. A fully expanded tree with the accession numbers for all genomes is shown in the Supplementary Material (Supplementary Figure S2). The scale bar indicates 10-amino-acid 10% sequence divergence. Highly supported nodes (SH-aLRT, • ≥ 95%, 1,000 replications) are indicated.
bathyarchaea BE326-BA-RLH (hereafter denoted as Bathy-6-B) and AD8-1 (hereafter denoted as Bathy-6-A). They are all highquality MAGs and were included in the subsequent analyses ( Table 1). Only bathyarchaeon SG8-32-3 (previously assigned to Bathy-1) was omitted because the completeness of the assembly (50.4%; based on our CheckM analysis) was too low for a reliable assessment of its metabolic capacity. Predicted genome sizes (1.0-2.5 Mbp), G++C contents (37.4-43.4 mol%), and coding densities (74.9-86.9%) of the MAGs from termite guts are in the same range as those of the other representatives of this subgroup ( Table 1). While the ANI values among the phylotypes of TB1 and TB2 range between 78.1 and 81.6%, the ANI values between members of TB1, TB2, and the other phylotypes of Bathy-6 are below the cutoff of the fastANI tool (<75%; Supplementary Figure S1), indicating that each lineage represents a separate genus-level taxon. This is confirmed by the results obtained with the Genome Taxonomy Database (GTDB) toolkit, which classified members of TB1 and TB2 as separate, genus-level lineages in the family UBA233 (order B26-1), a family that comprises also other members of Bathy-6. This indicates that TB1 and TB2 represent novel candidate genera in family UBA233, for which the names "Candidatus Termiticorpusculum" and "Candidatus Termitimicrobium" are proposed.
To identify the closest relatives of termite gut Bathyarchaeia and their respective habitats, we analyzed their phylogenetic position in the framework of rRNA genes available in public databases, which provides much better coverage than the small number of MAGs of the Bathy-6 subgroup available to date (Figure 2). The 16S rRNA gene sequences encoded by the MAGs form a well-supported monophyletic group with all other sequences of Bathyarchaeia that were previously obtained from the hindguts of higher termites (Friedrich et al., 2001;Shi et al., 2015;Grieco et al., 2019). Although each ribotype appears to be specific for a particular host species, the internal topology of the termite clade is not well resolved because of the large number of short sequences and the absence of 16S rRNA genes from many MAGs. The sequences in the termite clade are most closely related to clones obtained from a manure pit (EU662668; J. Ding, unpublished) and an anaerobic digestor fed with vinasses (U81774; Godon et al., 1997) and fall into the radiation of bathyarchaeal lineages in freshwater sediments, salt marshes, and anaerobic wastewater bioreactors (group 1.3 b; Ochsenreiter et al., 2003;Collins et al., 2005).

Capacity for CO 2 -Reductive Acetogenesis
We investigated the presence of all genes required for methanogenesis and reductive acetogenesis in all members of Bathy-6 with sufficiently complete genomes (Figure 3). All members of TB2 (phylotypes 8 and 9) encode the complete set of genes required for the reduction of CO 2 to acetyl-CoA via the archaeal version of the Wood-Ljungdahl pathway, using FIGURE 2 | 16S rRNA-based phylogeny of subgroup Bathy-6, indicating the placement of the termite clade among Bathyarchaeia from other environments. The maximum-likelihood tree is based on a curated alignment (1,424 positions) of all sequences in the SILVA database and their homologs retrieved from the bathyarchaeal MAGs and the low-quality bins obtained from the termite gut metagenomes (Hervé et al., 2020). The tree was rooted with members of Bathy-5 as outgroup. The scale bars indicate 0.05 nucleotide substitutions per site. SH-aLRT values (• ≥ 95%; • ≥ 80%, 1,000 replications) indicate node support. Branches marked with dashed lines indicate shorter sequences that were added using the parsimony tool. A fully expanded tree with the accession numbers of all sequences is shown in the Supplementary Material (Supplementary Figure S3). methanofuran (MFR) and tetrahydromethanopterin (H 4 MPT) as C 1 carriers (Figure 4). Formyl-MFR dehydrogenase is molybdenum-dependent (FmdABCDF; Hochheimer et al., 1996) and not the tungsten-dependent paralog. A homolog of fmdE, which occurs in methanogens, was not found in any of the MAGs, which suggests that the absence of subunit E is a characteristic feature of the bathyarchaeal complex. It has been shown that the Fmd complexes of Methanobacterium thermoautotrophicum and Methanosarcina barkeri are active also without this subunit (Hochheimer et al., 1996;Vorholt et al., 1996). Methylene-H 4 MPT dehydrogenase (Mtd) is more closely related to the NADH-dependent homolog of methylotrophic bacteria than to the F 420 H 2 -dependent homolog of methanogens. The CO dehydrogenase/acetyl-CoA synthase complex (CdhABCDE) and the (ADP-forming) acetyl-CoA synthetase (Acd; Musfeldt et al., 1999) are typical archaeal enzymes.
The same gene sets as in TB2 are also encoded by the more basal Bathy-6-S and Bathy-6-B (Figure 3), which indicates that the capacity to produce acetate from CO 2 might be a plesiomorphic trait of the Bathy-6 subgroup. The consistent absence of a key enzyme of the archaeal Wood-Ljungdahl pathway, methylene-H 4 MPT reductase (Mer), from all seven phylotypes (11 MAGs) of the TB1 lineage and from the most basal member of the subgroup, Bathy-6-A, suggests that the capacity to reduce CO 2 to the methyl level was lost at least twice during the evolutionary radiation of Bathy-6.
Homologs of the methyl-coenzyme M reductase (Mcr) complex, which encodes the key enzyme of methanogenesis, were not detected in any of the MAGs (Figure 4). Our observation contrasts with the report of Harris et al. (2018), who claimed that Bathy-6-B might represent an anaerobic methane oxidizer. However, their conclusion is based on the recovery of a 265-bp gene fragment classified as an mcrA gene in the original metagenome from which Bathy-6-B was assembled, i.e., not from the metagenomic bin. Considering also that the gene fragment FIGURE 3 | Gene functions encoded by termite gut bathyarchaea (TB1 and TB2) and other representatives of the Bathy-6 subgroup. All phylotypes with sufficiently complete genomes were included; their phylogenetic relationship was taken from Figure 1 (for strain designations, see Table 1). Colored circles indicate presence, and open circles indicate absence of the respective function; light blue indicates that a gene set is incomplete. The asterisks (*) in MtrABCDEFGH and FrhABG indicate that only MtrH or FrhB, respectively, is present. The number of aminotransferases encoded by each phylotype is indicated in the circle. If a phylotype is represented by more than one MAG, the annotation results were combined; details can be found in the Supplementary Material (Supplementary Table S2). H 4 MPT, tetrahydromethanopterin; MFR, methanofuran; Fpo, F 420 :methanophenazine oxidoreductase.
in question shows the highest similarity to a homolog from an uncultured euryarchaeal methanogen (GenBank: JX907770.1), it seems safe to conclude that members of the Bathy-6 subgroup are not methanogenic.
Although the capacity of Bathyarchaeia for reductive acetogenesis from CO 2 has been claimed repeatedly for several subgroups (He et al., 2016;Lazar et al., 2016;Yu T. et al., 2018;Zhou et al., 2018), the evidence was never fully conclusive. Actually, the comprehensive survey of all bathyarchaeal MAGs compiled by Zhou et al. (2018) lists only two MAGs that encode all genes required to operate the entire Wood-Ljungdahl pathway. One is the putatively methanogenic BA1 (Bathy-8) from a deep aquifer (Evans et al., 2015); the other is bathyarchaeon ex4484_135 (Bathy-15) from marine hydrothermal sediment (Dombrowski et al., 2017).

Capacity for Methylotrophic Acetogenesis
As all members of Bathy-6 encode a complete CO dehydrogenase/acetyl-CoA synthase (Cdh) complex (Figure 3), they might still synthesize acetyl-CoA using methyl groups derived from external sources. In all acetogenic bacteria and methylotrophic methanogens studied to date, the FIGURE 4 | Catabolic pathways encoded by the MAGs of termite gut Bathyarchaeia. The circles next to each enzyme indicate the presence of the coding genes in all, some, or none of the phylotypes of TB1 and TB2 (more details in Figure 3). Gray shading indicates pathways that are absent from all MAGs. The directionality of Fpo-like hydrogenase (Hfo) and ATP (synth)ase is discussed in the text. Dashed lines with question marks indicate hypothetical interactions. MFR, methanofuran; H 4 MPT, tetrahydromethanopterin (archaeal pathway); THF, tetrahydrofolate (bacterial pathway). A detailed list of genes present/absent in the respective MAGs is provided as Supplementary Material (Supplementary Table S2). methyltransferase systems consist of three components: (i) a set of substrate-specific methyltransferases (MT-I), (ii) their cognate methyl-accepting corrinoid proteins (CoP), and (iii) a second methyltransferase (MT-II) that transfers the methyl group of methyl-CoPs to THF (bacteria) or coenzyme M (archaea) (van der Meijden et al., 1983;Kreft and Schink, 1994;Kremp and Müller, 2020;Supplementary Figure S4A). We found that all MAGs of Bathy-6 encode CoPs that fall into the radiation of homologs assigned to other uncultured Archaea, with the CoPs of the di-and trimethylamine-specific methyltransferase systems (MtbC and MttC) of Methanomassiliicoccus luminyensis (Kröninger et al., 2017) and Acetobacterium woodii (Kremp et al., 2018) as closest relatives with a reliable functional annotation (Supplementary Figure S4). However, unlike the situation in methylotrophic bacteria and euryarchaea, where the CoP gene is colocalized with the gene of the cognate substrate-specific MT-I homologs (MtbB or MttB), the CoP gene of Bathy-6 is flanked by a gene encoding subunit H of tetrahydromethanopterin S-methyltransferase (MtrH; Supplementary Figure S4B).
In many methanogenic archaea, MtrH is part of the energy-conserving MtrABCDEFGH complex and catalyzes the transfer of the (CO 2 -derived) methyl group from methyltetrahydromethanopterin to the corrinoid prosthetic group of MtrA (Hippler and Thauer, 1999). However, in obligately methylreducing methanogens (Galagan et al., 2002;Borrel et al., 2014;Lang et al., 2015), which methylate CoM via their diverse methyltransferase systems (see above), the Mtr complex is absent. The presence of an isolated mtrH gene co-localized with a CoP gene has also been observed in the putatively methanogenic BA1 and BA2 (Bathyarchaeia) and several MAGs related to "Ca. Methanomethylicus mesodigestum" (Thermoproteota). It was proposed that the encoded proteins represent methyltransferase systems, which prompted the hypothesis that these uncultured lineages are methylotrophic methanogens (Evans et al., 2015;Vanwonterghem et al., 2016).
It is tempting to assume that also the CoP-MtrH couple of Bathy-6 is involved in the transfer of methyl groups from so far unidentified, substrate-specific methyltransferases to H 4 MPT (Figure 4). However, a catabolic role of the CoP-MtrH couple is not the only possible interpretation. In "Ca. Methanomethylicus mesodigestum, " the genes are colocalized with a homolog of metE encoding methionine synthase (Supplementary Figure S4B); it is also possible that the CoP-MtrH couple of Bathy-6 is involved in anabolic reactions that transfer methyl groups (provided by the cleavage of acetyl-CoA) from H 4 MPT to an unknown acceptor.

Hydrogen as Electron Donor
The operation of the Wood-Ljungdahl pathway requires electron donors in the form of reduced ferredoxin, NADH, and, in the case of archaea, also reduced cofactor F 420 (F 420 H 2 ) (Thauer et al., 2008;Schuchmann and Müller, 2014). The reduction of ferredoxin with H 2 is a critical step because it is endergonic at low hydrogen partial pressures and requires either an energyconverting hydrogenase or a flavin-based electron bifurcation system (Schut and Adams, 2009;Schuchmann and Müller, 2012).
Hydrogenases are present only in TB2 and the basal lineages of TB1 (Figure 3). One is a cytosolic, bidirectional [NiFe] hydrogenase of subgroup 3d, which uses NAD as electron acceptor . Phylogenetic analysis of the gene encoding the large subunit (hoxH) placed all homologs in a sister position to the Hox hydrogenases of phototrophic bacteria (Supplementary Figure S5). The gene order in the hoxEFUYH cluster is the same as in the gene clusters of other Hox complexes, which encode a prototypical heterodimeric [NiFe]-hydrogenase moiety (HoxHY) and a diaphorase moiety (HoxEFU); HoxEFU is homologous to the NuoEFG module of complex I and mediates the electron transport to NAD(P) (Eckert et al., 2012). Although members of group 3 are called "bidirectional hydrogenases, " hydrogen formation requires reduced ferredoxin or flavodoxin as electron donor (Gutekunst et al., 2014).
All MAGs that encode a Hox hydrogenase also possess a gene cluster that closely resembles those encoding the respiratory F 420 :methanophenazine oxidoreductases (Fpo) of Euryarchaeota and the homologous NADH:quinone oxidoreductases (Nuo/Nqo) of bacteria (complex I) ( Figure 5). As in other Fpolike or Nuo-like complexes, the genes encoding the FpoFO and NuoEFG modules, which provide substrate specificity for F 420 H 2 or NADH, respectively, are absent (Moparthi and Hägerhäll, 2011). However, six of the 11 subunits common to all Fpo and Nuo/Nqo complexes are also homologous to subunits of the energy-converting [NiFe] hydrogenases of group 4, underscoring their ancestral relationship to the respiratory complex I (Friedrich and Scheide, 2000;Schoelmerich and Müller, 2019).
Classification with HydDB placed the D subunit of the 11subunit complex of the Bathy-6 MAGs among the catalytic subunits of [NiFe] hydrogenases in subgroup 4 g. The hydrogenases in subgroup 4 g are structurally heterogeneous and differ fundamentally both in the number of their subunits and the arrangement of their coding genes Schoelmerich and Müller, 2019; Figure 5). Their large subunits form several distinct phylogenetic lineages (Subgroups 4g-1 to 4g-6; Figure 6), which indicates that they evolved independently from each other. The gene cluster encoding the Fpo-like hydrogenase complex of Bathy-6 (hereafter referred to as Hfo) has an organization almost identical to that of the corresponding clusters of Ca. Methanomethylicus mesodigestum (Thermoproteota) and Pyrodictium delaney (Crenarchaeota) (Figure 5), whose large subunits represent phylogenetic sister groups of subgroup 4g-6 ( Figure 6). The coordination sites of the [NiFe] cofactor on the large subunit of all [NiFe] hydrogenases (L1 and L2 motifs; Vignais and Billoud, 2007), which are no longer conserved in NuoD and FpoD, are present in all Bathy-6 homologs (Figure 7).
The Hfo hydrogenase of Bathyarchaeia is most interesting from an evolutionary perspective, as it represents the first [NiFe] hydrogenase that is composed of the same 11 subunits and shares the same organization of the coding genes as the archaeal Fpo complex and the bacterial Nuo/Nqo complex (Figure 5). Both Hfo and the predicted Hfo-like complexes of Pyrodictium delaneyi and Ca. Methanomethylicus mesodigestum (subgroup 4g-6) lack the Na + transport module of the membrane-bound hydrogenase (Mbh) complex of Pyrococcus furiosus (subgroup 4d), whose similarity to the respiratory complex I (Nqo) of Thermus thermophilus has been well documented (Yu H. et al., 2018). Notably, the Na + transport module (MbhABCF) in the Mbh of P. furiosus is also present in the gene cluster encoding the [NiFe] hydrogenase of Thermosphaera aggregans (subgroup 4g-5) and other members of Desulfurococcales (not shown), which encode all subunits of the Mbh complex of P. furiosus, albeit in a different gene order. The striking synteny between the gene clusters encoding the Hfo of Bathyarchaeia and the Fpo-like complex of Methanomassiliicoccales, including the absence of genes encoding the F 420 -binding module (FpoFO), and the phylogeny of its large subunit suggest that the Hfo complex represent a closer evolutionary link between the energyconverting hydrogenases and the modern respiratory complexes than the Mbh of Thermococci.
None of the hydrogenases of subgroup 4 g have been biochemically characterized, but they are presumed to couple FIGURE 5 | Organization of the gene clusters encoding the respiratory complexes (blue background) and the ancestral [NiFe] hydrogenases in group 4 (red background). Identical colors indicate homologous genes; a phylogenetic analysis of the catalytic subunit of [NiFe] hydrogenases and its homologs (red) is shown in Figure 6. The font style of the gene labels indicates differences in the subunit nomenclature of Nuo/Fpo (uppercase), Mbh (lowercase), and Ech (italics).
Frontiers in Microbiology | www.frontiersin.org • ≥ 80%, 1,000 replications) indicate node support. The genomic context of the highlighted genes is shown in Figure 5. Gene numbers indicate IMG/Mer gene IDs.
the formation of H 2 from reduced ferredoxin to the formation of an electrochemical membrane potential Søndergaard et al., 2016;Schoelmerich and Müller, 2019). This is in agreement with biochemical data obtained for the Fpo-like 11-subunit complex of methanogenic Euryarchaeota, which generate an electrochemical membrane potential during electron transport from reduced ferredoxin to methanophenazine (Methanosaeta; Welte and Deppenmeier, 2011) or a so far unidentified electron acceptor (Methanomassiliicoccales; Kröninger et al., 2016). The absence of genes involved in the biosynthesis of methanophenazine from all MAGs of Bathy-6 (Supplementary Table S2) adds to the evidence that the Hfo of Bathyarchaeia is not a respiratory complex but is instead a novel energy-converting hydrogenase that catalyzes the reduction of ferredoxin with H 2 using the electrochemical membrane potential (Figure 4).
While Hox and Hfo hydrogenase should provide members of TB2 with the NADH and reduced ferredoxin required to operate the Wood-Ljungdahl pathway, the source of F 420 H 2 as potential electron donor for methylene-H 4 MPT reductase (Mer) remains unclear. All phylotypes encode enzymes involved in the biosynthesis of F 420 (Supplementary Table S2), but a complete gene set encoding F 420 -reducing [NiFe] hydrogenase (FrhABG, subgroup 3a; Supplementary Figure S5) is present only in Bathy-6-A. All members of TB2 and several phylotypes of TB1 encode a homolog of FrhB, an iron-sulfur flavoprotein with an F 420 -binding site, but not the hydrogenase subunits (Figure 3). It is possible that FrhB is involved in the reduction of F 420 via an interaction with HdrABC and an unknown electron donor, as proposed for the methane-oxidizing Ca. Methanoperedens spp. (Arshad et al., 2015). The only member of subgroup Bathy-6 that encodes a complete FrhABG is Bathy-6-A. It is also the only MAG that encodes a methylviologen-dependent [NiFe] hydrogenase (MvhADG; Supplementary Figure S5, subgroup 3c), which forms an electron-bifurcating complex with the soluble heterodisulfide reductase (HdrABC) and catalyzes the hydrogendependent reduction of ferredoxin and the heterodisulfide of coenzyme M (CoM) and coenzyme B (CoB) in methanogens (Kaster et al., 2011). The presence of genes encoding HdrABC, MvhADG, and a complete Wood-Ljungdahl pathway in the putatively methanogenic BA1 (Bathy-3) provides strong evidence that BA1 is capable of hydrogenotrophic methanogenesis (Evans et al., 2015). In Bathy-6-A, however, the pathway is incomplete, and the identity of the heterodisulfide reduced by Hdr remains unclear. Interestingly, the same constellation as in Bathy-6 has been recently reported for the bathyarchaeal MAG CR_14 from marine sediments, which represents another, novel subgroup of Bathyarchaeia (Farag et al., 2020).

Organic Substances as Electron Donors
Most members of TB1 and all basal lineages of Bathy-6 lack Hox and Hfo (Figure 3), which means that they cannot grow lithotrophically with H 2 as electron donor. However, the reduced Fd required to operate reductive acetogenesis, either via the Wood-Ljungdahl pathway (TB2) or by methylotrophy (all phylotypes), could be provided also by the oxidation of organic substrates (Figure 4). Such organotrophic acetogenesis is common among bacteria with a homoacetogenic lifestyle (Drake, 1994;Schink, 1994). All Bathy-6 genomes (except Bathy-6-A) encode pyruvate:ferredoxin oxidoreductase (Por) and indolepyruvate:ferredoxin oxidoreductase (Ior), and some also encode 2-oxoglutarate:ferredoxin oxidoreductase (Oor), all of which catalyze the oxidative decarboxylation of 2-oxo acids to their corresponding acyl-CoA esters (Figure 3). The 2-oxoacids would result from the transamination of amino acids via numerous aminotransferases encoded by all genomes; a putative amino acid permease, however, is encoded only in TB2. ATP would be formed via the ADP-dependent acetyl-CoA synthetase, which accepts also other acyl substrates in P. furiosus . Such pathways have been shown to operate in other archaea (P. furiosus, Thermococcus spp.; Kengen and Stams, 1993;Heider et al., 1996) and in the insect gut-associated bacterium Elusimicrobium minutum (Herlemann et al., 2009) during growth on glucose, where they result in a net formation of alanine.
The data compiled by Zhou et al. (2018) suggest that several lineages of Bathyarchaeia, including Bathy-6-A; Lazar et al., 2016), have the capacity to ferment various organic carbon compounds. However, genes encoding extracellular peptidases, which are numerous in other Bathyarchaeia, seem to be less prevalent in the MAGs of Bathy-6 and Bathy-1 (Feng et al., 2019), which suggests that members of these subgroups are limited to the utilization of amino acids or oligopeptides that are small enough to be transported across the cytoplasmic membrane.
There is no indication that members of Bathy-6 have the capacity to utilize sugars. Like Bathy-6-A (Lazar et al., 2016), all MAGs of TB1 and TB2 encode many genes of the classical Embden-Meyerhof-Parnas (EMP) pathway, including glyceraldehyde-3-phosphate dehydrogenase and phosphoglycerate kinase. However, all MAGs lack hexokinase and the alternative archaeal glycolytic enzymes (Bräsen et al., 2014), and most MAGs lack phosphofructokinase and pyruvate kinase. As all MAGs encode phosphoenolpyruvate synthetase and fructose bisphosphatase, it is likely that the EMP pathway functions only in gluconeogenesis. Sugar transporters were not detected; the role of the lipooligosaccharide ABC transporter encoded by almost all phylotypes from termite guts (except phylotype 9) is not clear (Supplementary Table S2). The identification of a cellulolytic system in Bathy-6-A (Lazar et al., 2016) requires verification.

Energy Conservation in TB2
In acetogenic bacteria growing on hydrogen and CO 2 , all ATP synthesized by SLP is consumed in the activation of formate. Therefore, energy conservation involves electron-transport phosphorylation, which is driven by the oxidation of reduced ferredoxin via membrane-bound electron-transport complexes (Schuchmann and Müller, 2014;Basen and Müller, 2017). By contrast, the activation of formate (i.e., the formation of formylmethanofuran) in the archaeal variant of the Wood-Ljungdahl pathway is not ATP-dependent but is instead driven by the reducing power of ferredoxin, yielding a full ATP per acetate produced via SLP. However, thermodynamics dictates that a fraction of this ATP must be reinvested, as a metabolism where the net ATP yield exceeds the free-energy change of the reaction would become endergonic (Thauer et al., 2008).
Fermenting bacteria that lack respiratory chains energize their membrane by operating their ATP synthase in the reverse direction (Buckel and Thauer, 2013). Likewise, members of Bathy-6 that possess a complete Wood-Ljungdahl pathway (i.e., the phylotypes in TB2) might use part of the ATP gained by SLP to generate an electrochemical membrane potential that drives the H 2 -dependent reduction of ferredoxin via Hfo (see above). Other energy-converting complexes that would allow generation of reduced ferredoxin, such as the group-4 [NiFe] hydrogenases in acetogenic bacteria and methanogenic archaea (Ech, Künkel et al., 1998;Eha and Ehb, Tersteegen and Hedderich, 1999) or an NADH:Fd oxidoreductase complex (RnfABCDEG, Westphal et al., 2018), were not detected in any member of Bathy-6. If one assumes that the ATPase translocates 4 H + per ATP and Hfo translocates only 2 H + during electron transport from H 2 to ferredoxin (m = 4, n = 2 in Figure 4), production of 2 Fd 2− via Hfo would completely consume the energy conserved by SLP (1 ATP). Therefore, it is likely that members of TB2 grow mixotrophically, producing one Fd 2− from H 2 (via Hfo) and the other by the oxidation of pyruvate or other 2-oxo acids. An entirely lithotrophic pathway would be feasible if one Fd 2− is produced by Hfo and the other, together with F 420 H 2 , by flavin-based electron bifurcation (see above), but this would require an additional, unknown electron donor.
It is intriguing that several phylotypes of TB1 and TB2 (Figure 3) and also bathyarchaeal MAGs from other subgroups (Evans et al., 2015;Zhou et al., 2018) do not encode an ATP synthase (neither the genes for the archaeal V-type ATP synthase nor those for the bacterial equivalent were detected). While this observation is most likely explained by incomplete genome assemblies, it cannot be entirely excluded that these organisms generate their membrane potential (vital for any organism) by other means. In this case, the Hfo complex (if present) might operate in the reverse direction, using reduced ferredoxin provided by the oxidation of organic substrates to produce H 2 and generate an electrochemical membrane potential, like the energy-converting hydrogenases in fermenting bacteria.
In principle, the entire Wood-Ljungdahl pathway is reversible and can oxidize acetate to CO 2 given the appropriate thermodynamic framework. This has been demonstrated in syntrophic cultures of "Reversibacter"-like microorganisms with hydrogenotrophic partners (Lee and Zinder, 1988;Schnürer et al., 1997) and has been suggested to occur also in Bathyarchaeia (Evans et al., 2015;Xiang et al., 2017). However, at least in the termite hindgut, where the hydrogen partial pressure is much higher than in sediments (Ebert and Brune, 1997;Schmitt-Wagner and Brune, 1999) and reductive acetogenesis often prevails over methanogenesis as electron sink (Brauman et al., 1992;Tholen and Brune, 1999;Tholen and Brune, 2000), an anaerobic oxidation of acetate is an unlikely scenario.

Ecological Aspects
Although the proportion of archaeal rRNA in termite hindguts is relatively small (0.9-2.3% of all prokaryotic rRNA; Brauman et al., 2001), methanogenesis represents a substantial hydrogen sink (Brune, 2019). Considering that the proportion of reads assigned to bathyarchaeal MAGs in the hindgut metagenomes of higher termites (0.03-2.5%; avg. 0.69%) is four times higher than that assigned to euryarchaeal MAGs (0.02-0.79%; average, 0.16%; Supplementary Table S2 in Hervé et al., 2020), the population sizes of Bathyarchaeia might be sufficient to contribute significantly to acetogenesis, particularly in soilfeeding species.
However, the substrates of termite gut Bathyarchaeia remain open to speculation. While only members of TB2 have the genomic capacity for lithotrophic acetogenesis, almost all members of Bathy-6 have the capacity to ferment amino acids and might employ organotrophic acetogenesis from methylated substrates as an electron sink. This would explain their prevalence in soil-and humus-feeding termites. It has been estimated that soil peptides and other nitrogen-rich humus constituents contribute substantially (20-40%) to the dietary carbon oxidized by soil-feeding Cubitermes spp. (Ngugi et al., 2011), which is consistent with the depletion of peptides in soil organic matter during gut transit (Griffiths et al., 2012) and the high ammonia concentrations (up to 130 mM) in the posterior hindgut (Ji and Brune, 2006). An NifDH homolog (pfam00142 and pfam001428) encoded by both TB1 and TB2 is most likely not involved in dinitrogen fixation but rather in a so far unidentified archaeal tetrapyrrole biosynthesis pathway (Ghebreamlak and Mansoorabadi, 2020).
Stable-isotope probing of salt marsh sediments indicated that members of Bathy-8 and Bathy-6 assimilate organic substrates, notably excluding proteins and inorganic carbon (Seyler et al., 2014). Yu T. et al. (2018), however, reported that the addition of lignin to an estuarine sediment sample selectively stimulated the growth of Bathy-8 and the incorporation of carbon from 13 C-bicarbonate into archaeal tetraether lipids, which suggests that members of Bathy-8 are methylotrophs that use ligninderived methyl groups. Together with the potential capacity for methyl group utilization in many bathyarchaeotal MAGs (Seyler et al., 2014;Yu T. et al., 2018; this study), these results explain the observations of Lever et al. (2010), who found that porewater acetate in deep-subseafloor sediments was depleted in 13 C relative to sedimentary organic matter and postulated that a substantial fraction of the acetate produced in marine sediments might stem from reductive acetogenesis, fueled by microbial fermentation products, molecular hydrogen, and the methoxy groups of lignin monomers.
The utilization of the methoxy groups of lignin-derived aromatic compounds is a common trait of many acetogenic bacteria (Schink et al., 1992;Drake, 1994). Methoxylated aromatic compounds are demethylated by the hindgut microbiota of termites (Brune et al., 1995), but the organisms responsible for this activity have not been identified. It is tempting to speculate that termite gut Bathyarchaeia are organotrophic (TB1) or mixotrophic (TB2) acetogens that utilize methylated compounds such as lignin derivatives as methyl group donors and reduce CO 2 either with molecular hydrogen and/or with reducing equivalents derived from the oxidation of organic substrates.
It has been speculated that acetogenic archaea might have an energetic advantage over acetogenic bacteria, as they do not have to invest ATP to activate formate (He et al., 2016). However, the net synthesis of ATP is limited by the free-energy change of an acetogenic metabolism, which is independent of its reaction path and requires part of the ATP gained by SLP to be reinvested (e.g., for ferredoxin reduction; see above). Rather, it is feasible that the capacity for methylotrophic acetogenesis, which is less sensitive to low hydrogen partial pressures than hydrogenotrophic acetogenesis, provides an energetic advantage, analogous to the situation in methyl-reducing methanogens (Feldewert et al., 2020). Moreover, it has been argued that long generation times contribute to the difficulties surrounding the enrichment and isolation of Bathyarchaeia in the laboratory (Yu H. et al., 2018). In view of the relatively short residence time of organic matter in termite guts (24-48 h;Kovoor, 1967;Bignell et al., 1980), the growth rates of termite gut Bathyarchaeia must be high enough to avoid washout -unless they are attached to the intestinal surface.
Habitat: The hindgut of higher termites.
Habitat: The hindgut of higher termites.

CONCLUSION
To date, the nonmethanogenic archaea in termite guts and their potential role in symbiotic digestion have received little attention. Our study provides strong evidence that termite gut Bathyarchaeia and other members of the Bathy-6 subgroup are archaeal acetogens; they possess the genomic potential to conserve energy by the production of acetyl-CoA from CO 2 (Ca. Termitimicrobium; TB2) and/or possibly methyl groups (almost all members of Bathy-6, including Ca. Termiticorpusculum; TB1). As in bacterial acetogens, their energy metabolism is likely mixotrophic or organotrophic. We identified a complete gene set encoding a novel Fpolike 11-subunit hydrogenase, which closes the evolutionary gap between the ancestral [NiFe] hydrogenases and the respiratory complex I and would enable members of TB2 to grow mixotrophically on H 2 . All members of Bathy-6 are probably able to derive reducing equivalents from the oxidation of organic substrates (viz., amino acids) and use reductive acetogenesis as an electron sink. These findings agree with previous claims concerning the capacity for reductive acetogenesis in other subgroups of Bathyarchaeia. However, this is the first time that all genes encoding the Wood-Ljungdahl pathway and the components required for the provision of reducing equivalents and energy conservation are conclusively documented. Although eight of the nine closely related phylotypes of termite gut Bathyarchaeia were represented by high-quality MAGs, a complete pathway was detected only in members of TB2 and two more basal lineages from other environments. This underscores the longstanding caution that the mere presence of marker genes of the Wood-Ljungdahl pathway does not qualify an organism as an acetogen, as many of its enzymes are found also in nonacetogenic organisms, where they are involved in the assimilation and interconversion of C 1 metabolites (Drake, 1994).

Metagenome-Assembled Genomes
Data on the MAGs from termite guts are from Hervé et al. (2020). All other MAGs were retrieved from the NCBI Assembly database 1 ; accession numbers are listed in Table 1. Assembly coverage was determined as described by Hervé et al. (2020). Average nucleotide acid identities (ANIs) were calculated with fastANI (Jain et al., 2018). Protein-coding genes were predicted with Prodigal v2.6.3 (Hyatt et al., 2010).

16S rRNA Gene Phylogeny
SSU rRNA gene sequences in the MAGs and other bathyarchaeotal bins obtained from the original metagenomes (Hervé et al., 2020) were identified using the ssu_finder function implemented in CheckM. Sequences were imported into the alignment of rRNA gene sequences in the SILVA SSURef NR database release 132 ( 3 Quast et al., 2013) using Arb v6.0.6 (Ludwig et al., 2004). After automatic alignment of the imported sequences using the PT server and the Fast Aligner tool implemented in Arb, the alignment was manually refined using the Arb editor, considering secondary structure information to identify homologous base positions. After removing sites with more than 50% gaps, the alignment consisted of 1,424 sites with unambiguously aligned base positions. Phylogenetic trees were reconstructed by maximum-likelihood analysis with IQ-TREE using the best-fit evolutionary model (GTR+F+R4) suggested by ModelFinder; node support was assessed using SH-aLRT with 1,000 resamplings. Gene fragments (<1,300 bp) were inserted into the core tree using the parsimony tool implemented in Arb.

Gene Discovery and Annotation
For an initial exploration of the genes potentially involved in energy metabolism, bathyarchaeotal MAGs were analyzed using the annotation provided in the IMG/Mer database ( 4 Chen et al., 2019). Annotation results were verified, and missing functions were identified with hidden Markov model (HMM) searches, using HMMER v3.1b2 (Eddy, 2011) with a threshold E-value of 1E-5; the respective models are listed in Supplementary  Table S3. The identity of all genes of interest was confirmed using the NCBI Conserved Domain search (Marchler-Bauer and Bryant, 2004) and BLASTp (Altschul et al., 1990). Additionally, Bathy-6-S and Bathy-6-B were annotated with BlastKOALA (Kanehisa et al., 2016). When indicated, closest neighbors were identified by BLAST and aligned using MAFFT v7.305b with the L-INS-i method (Katoh and Standley, 2013). Phylogenetic trees were reconstructed by maximum-likelihood analysis with IQ-TREE (Nguyen et al., 2015) using the best-fit evolutionary model (LG+G+I) suggested by ModelFinder (Kalyaanamoorthy et al., 2017). Node support was assessed using SH-aLRT with 1,000 resamplings (Anisimova et al., 2011).

Analysis of [NiFe] Hydrogenases
Putative [NiFe] hydrogenase genes were identified by HMM searches (see above), using the highly resolved models provided by Anantharaman et al. (2016). Search results were confirmed with HydDB, a web-based tool for hydrogenase classification and analysis ( 5 Søndergaard et al., 2016).
The deduced amino acid sequences of the large subunit (LSU) of [NiFe] hydrogenases recovered from the MAGs and their top BLAST hits on the IMG/Mer database were imported into an alignment of NuoD and FpoD homologs (Lang et al., 2015), which was completed with representative members of other hydrogenase classes extracted from HydDB. The alignment was manually refined in the Arb editor. Phylogenetic trees were reconstructed by maximum-likelihood analysis with IQ-TREE (Nguyen et al., 2015) using the best-fit evolutionary model (LG+G+I) suggested by ModelFinder (Kalyaanamoorthy et al., 2017). Node support was assessed using SH-aLRT with 1,000 resamplings (Anisimova et al., 2011).

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
HL and AB designed the study. HL analyzed data and wrote the first draft of the manuscript. VH contributed to the analyses. AB analyzed data and revised the manuscript. All authors edited and approved the final version of the manuscript.

ACKNOWLEDGMENTS
The authors thank the Joint Genome Institute for their metagenome sequencing service and for providing the IMG/ER platform.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020. 635786/full#supplementary-material Supplementary Figure 1 | Average nucleotide identity (ANI) of the MAGs in subgroup Bathy-6. The termite gut Bathyarchaeia were assigned to phylotypes based on ANI > 99%. NA indicates ANI values <75%, which are not returned by the fastANI program.
Supplementary Figure 2 | Genome-based phylogeny of termite gut Bathyarchaeia illustrating the relationship of lineages TB1 and TB2 to other MAGs in the Bathy-6 subgroup. MAGs mentioned in the text are marked in bold. The maximum-likelihood tree was inferred from a concatenated alignment of 43 proteins using the LG+F+I+G4 model and rooted with selected Crenarchaeota and Euryarchaeota as outgroup. The numbers in circles indicate the phylotypes discussed in the text ( Table 1). MAGs included in the comparative analysis ( Figure 3) are shown in bold. The tree was rooted other archaeal genomes as outgroup. The scale bar indicates 10-amino-acid substitutions per site. Node support values (SH-aLRT) are shown in blue. A simplified version of the tree is shown in Figure 1.
Supplementary Figure 3 | 16S rRNA-based phylogeny of subgroup Bathy-6, indicating the placement of the sequences from termite guts among those obtained from other environments. The maximum-likelihood tree is based on a curated alignment (1,424 positions) of all sequences in the SILVA database and their homologs retrieved from the bathyarchaeal MAGs (in bold) and the low-quality bins obtained from the termite gut metagenomes (Hervé et al., 2020). The tree was rooted using members of Bathy-5 as outgroup. The scale bars indicate 0.05 nucleotide substitutions per site. Node support values (SH-aLRT) are shown in blue. Branches marked with dashed lines indicate shorter sequences that were added using the ARB parsimony tool. A simplified version of the tree is shown in Figure 2. Supplementary Table 1 | Taxonomic assignment and characteristics of the bathyarchaeotal MAGs from termite guts (from Hervé et al., 2020).
Supplementary Table 2 | Annotation details of the genes that encode the metabolic pathways and other functional markers in the 15 bathyarchaeotal MAGs from termite guts, as discussed in the text (see Figures 3, 4).