Integrative modeling of gene and genome evolution roots the archaeal tree of life

Identifying a Consensus Unrooted Topology for the Archaea with Concatenated Proteins and Multigene Supertrees.

We used the OMA algorithm (22) to identify single-copy orthologs on 62 archaeal genomes sampled from across the known diversity of the domain. Our sample included 21 single-cell genomes and metagenomic bins from uncultivated lineages, which are now known to represent some of the most abundant and ecologically important Archaea (4, 5). We filtered candidate marker genes to remove potential horizontal gene transfers (HGTs) (Materials and Methods) and inferred a concatenated protein phylogeny (Fig. 1) for a Dayhoff-recoded (23) supermatrix comprising 45 proteins under the CAT+GTR (generalized time-reversible) model, the best-fitting evolutionary model. To complement the supermatrix analysis, we used a multigene supertree approach, matrix representation with parsimony (MRP) (24–27), which integrates the phylogenetic signal for vertical descent from a much broader sample of genes than can be accommodated by concatenation alone (Materials and Methods). The supertree inferred by MRP fitted to a dataset of 3,242 single gene trees (SI Appendix, Fig. S1) is in good agreement with the concatenation tree topology (Fig. 1), providing a robust phylogenetic backbone for rooting analysis and suggesting strong vertical signals in the data. The main difference between the two trees is the position of the Thermococcales, which emerge at the base of a clade comprising the TACK Archaea and Lokiarchaeum in the concatenation tree but at the base of the Euryarchaeota in the supertree. This clade has been difficult to place in previous analyses, and both of the positions that we recovered have been supported by previous work (9, 10, 18, 28). We evaluated both positions for the Thermococcales in our rooting analysis.

The unrooted phylogeny contains three major clans (29), or potential clades defined by a single split on the tree. These correspond to (i) a metabolically diverse assemblage comprising the TACK (8) Archaea—Thaumarchaeota (2), Aigarchaeota (30), Crenarchaeota (31), Korarchaeota (32)—and the recently discovered Lokiarchaeota (6), which emerges as the sister group of TACK; (ii) the core Euryarchaeota, comprising the methanogenic Euryarchaeota and their relatives (31), with the exception of the Thermococcales; and (iii) the DPANN Archaea (4), comprising Nanoarchaeum and its relatives. Many of the DPANN Archaea described to date have small cells (<1 μm) (5) and reduced genomes in which many core metabolic pathways are incomplete (5). The first DPANN to be characterized, Nanoarchaeum equitans, is obligately dependent on the crenarchaeote Ignicoccus hospitalis for growth (33), and other members of the group also have been observed in direct contact with larger archaeal cells (34), suggesting that symbiotic or parasitic lifestyles may be a common feature of the DPANN lineage. Nevertheless, genome analyses suggest that at least some members of the group might be capable of a free-living lifestyle (35).

Are the DPANN Archaea a Clan?

Our phylogeny is in agreement with recent published analyses in recovering the clanhood of TACK (8, 28, 36, 37) and the core Euryarchaeota (5, 9, 10), and supports the placement of Lokiarchaeum at the base of the TACK phylum (6). DPANN clanhood is in agreement with some recent reports (4, 5, 38), but has been challenged (9) on the grounds that the high rates of sequence evolution shared by some DPANN lineages make them vulnerable to long-branch attraction (LBA) (39, 40). DPANN clanhood raises unsettling parallels with early molecular phylogenies of the eukaryotes, in which fast-evolving parasitic lineages were drawn to the base of the tree by LBA (41). Published DPANN phylogenies have shown conflicting results (2, 5, 9, 42), although analyses using the CAT+GTR model, which may be less susceptible to LBA than simpler methods (43), have recovered DPANN monophyly (38). Owing to this uncertainty, recent analyses of the archaeal root have excluded DPANN (10) or have included only a subset of sequenced lineages (9), on the grounds that their presence would interfere with the overall resolution of the tree. The limitation of this approach is that DPANN lineages are ecologically important (44) and represent a substantial proportion of known archaeal diversity (4, 5). Therefore, any analysis of the archaeal root that does not account for their origins is necessarily incomplete.

We performed a series of tests designed to investigate whether DPANN clanhood could be attributed to LBA. We first recoded the alignment into the four Dayhoff categories (23, 45), which made the data easier to model by reducing both compositional heterogeneity and substitutional saturation (46). Analysis of this recoded matrix under CAT+GTR, one of the best phylogenomic models for ameliorating the effects of LBA (43), recovered DPANN clanhood with maximal posterior support [posterior probability (PP) = 1; Fig. 1]. We then selectively removed the longest-branching DPANN lineages from the analysis (Materials and Methods and SI Appendix, Fig. S2), and again obtained maximal support for DPANN clanhood, although Pacearchaeota now clustered within Woesearchaeota with moderate support (PP = 0.89). We also reanalyzed the original alignment after applying a more stringent approach to identify and remove fast-evolving sites (the BLOSUM62 matrix in BMGE) (47), which are considered the sites most susceptible to LBA artifacts (16); support for DPANN clanhood was unchanged (SI Appendix, Fig. S3).

Next, we reasoned that if the DPANN lineages were being artifactually drawn to the base of the tree because of LBA, then an analysis of DPANN alone might not reveal the same in-group topology as that seen in the full analysis, including the euryarchaeotal and TACK outgroups (17). Artifacts of this type have previously been observed in analyses of within-eukaryote relationships, whereby fast-evolving eukaryotes that were drawn toward the prokaryotic outgroup were recovered in the expected position in a eukaryote-only analysis (fig. S62 in ref. 48). However, a CAT+GTR analysis of the DPANN portion of the concatenation alone resulted in a topology compatible with that of the overall tree (SI Appendix, Fig. S4), with the exception of the position of Haloredivivus, which exchanges with its nearest neighbor in the DPANN-only reanalysis at PP = 0.88. We also considered the possibility that DPANN clanhood might be an artifact of nonrandom gene representation in the supermatrix, which potentially could lead to systematic error (49). Because many DPANN lineages were represented by an incomplete metagenomic bin, they often contained more gene absences (and thus gaps) in the supermatrix compared with other Archaea (DPANN gene representation ranging from 4 to 40 genes; median, 28 genes). To evaluate the impact of this gene representation bias on our analyses, we subsampled the original dataset, selecting the most complete DPANN genomes (10 genomes, including at least one genome from each major DPANN sublineage) and the most widely conserved genes (25 genes) to produce a supermatrix in which gene representation was equal across Euryarchaeota, TACK, and DPANN. Analysis of this supermatrix under CAT+GTR resulted in a topology (SI Appendix, Fig. S5) almost identical to that inferred from the full concatenation (Fig. 1). In particular, support for both the in-group relationships within DPANN and the clanhood of the group as a whole were identical to those seen in the original analysis.

Finally, and because DPANN Archaea tend to have above-average evolutionary rates, we considered the possibility that the apparent clanhood of the group as a whole is the result of LBA between the stems leading to each distinct sublineage. Thus, we investigated the behavior of individual sublineages in a series of concatenated protein analyses from which all other DPANN Archaea had been removed (SI Appendix, Figs. S6–S12). The idea is that if DPANN are monophyletic, then their constituent lineages should each individually connect to the same point on a subtree containing only Euryarchaeota and the TACK/Lokiarchaeum clade. In these analyses, the Diapherotrites, Aenigmaarchaeota, and Woesearchaeota lineages—composing just over one-half (14 of 24) of sampled DPANN lineages—fell between the Euryarchaeota and TACK/Lokiarchaeum clans, as in the full analysis; however, the remaining DPANNs grouped at the base of the Euryarchaeota, either with the Thermococcales (Nanoarchaeum) or within the cluster 1 methanogens (Nanohaloarchaeota, Pacearchaeota, and the solitary Parvarchaeum). The difficulties in finding stable positions for single DPANN lineages including Nanoarchaeum, Nanohaloarchaeota, and Parvarchaeum is already clear from comparing trees in previously published work (2, 10, 38, 42, 50). The results of single-lineage analyses are also difficult to compare with the full analysis, because there is no principled statistical framework within which to evaluate whether attempting to place DPANN lineages individually ameliorates or aggravates potential phylogenetic artifacts, such as LBA. Better taxonomic sampling has been shown to improve phylogenetic inference (51–53), and there is no posterior support for these alternative placements from any of our analyses in which DPANN monophyly was tested directly, including our supertree analysis and the series of supermatrix analyses performed with methods commonly used to ameliorate LBA (23, 43, 54) (SI Appendix, Figs. S2–S5). Nonetheless, we considered both possibilities—monophyletic and polyphyletic DPANN—and also performed an analysis in which all DPANN lineages were excluded in our subsequent rooting and gene content analyses.

Using a Bacterial Outgroup to Root the Archaea.

Recent work using a bacterial outgroup to root the Archaea has recovered a root either between Euryarchaeota and TACK (9) or within the Euryarchaeota (10). Our own outgroup rooting analysis using 29 universally conserved protein-coding genes and the CAT+GTR model (SI Appendix, Table S2 and Fig. S13) did not robustly distinguish between these two hypotheses. We obtained weak to moderate posterior support for the exclusion of four clades from the root (DPANN, TACK/Lokiarchaeum+Thermococcales, core Euryarchaeota, and cluster I methanogens), whereas the basal split within the Archaea was unresolved. The outgroup approach allows the addition of a priori rooting information to trees inferred under standard models of sequence evolution, which do not directly infer the root (55, 56). However, outgroup rooting is known to be problematic when the outgroup is distantly related to the in-group (16, 17, 57), as is the case when one cellular domain is used to root another. The length of the branch leading to the outgroup is particularly striking in our analysis (SI Appendix, Fig. S14), where the bacterial stem was predicted to have experienced 4.79 substitutions per site, compared with a mean of 0.192 (range, 0.0157–0.546) for within-domain branches. The use of long outgroup branches is a general problem that has contributed to disagreements about the archaeal root as well as about the roots of other major radiations (58–61), motivating a search for alternative rooting methods.

Bringing More Data to Bear on the Archaeal Root.

The use of gene duplications to root major clades has a venerable history in molecular evolution (45), particularly for resolving the root of the tree of life (11–13, 62–64). More recently, it has been appreciated that gene gains, losses, and horizontal transfers also contain information about the root of a species tree that can be integrated using probabilistic gene tree-species tree reconciliation approaches (65–67). We used a recently developed method known as amalgamated likelihood estimation (ALE) (67) to calculate gene family likelihoods for each of the 31,236 homologous gene families encoded by our sample of 62 archaeal genomes, under a set of candidate root positions on the archaeal species tree (Fig. 1) corresponding to published rooting hypotheses as well as a selection of other plausible rooting positions, such as between each of the major lineages (SI Appendix, Tables S3–S5). We also evaluated a species tree in which the DPANN were polyphyletic (SI Appendix, Fig. S15), as has been suggested by some single-lineage supermatrix analyses. Different roots on the species tree imply different scenarios of gain, duplication, transfer, and loss for the gene families observed on modern genomes (Fig. 2), and because of this they have different likelihoods under the model. The rates of gene duplication, transfer, and loss (DTL) were inferred from the data using maximum likelihood (ML) optimization, and ALE incorporates uncertainty in the underlying gene trees using conditional clade probabilities (68, 69). This means that poorly supported disagreements between the species tree and the gene trees do not unduly affect estimates of the number of DTL events.

Using an approximately unbiased test (70) to establish a confidence set from our analyses for the archaeal root at P > 0.05, we were able to reject all but a single root (SI Appendix, Table S3), the root between DPANN and a clade comprising the Euryarchaeota and TACK/Lokiarchaeum lineages (Fig. 3A). We obtained significantly higher likelihoods using rooted trees in which the Thermococcales were placed at the base of the Euryarchaeota rather than with the TACK Archaea, in agreement with our supertree (SI Appendix, Fig. S1) and some previous analyses (9, 10, 48). The tree in which DPANN were polyphyletic had the worst likelihood score of the trees considered and was rejected at P = 4 × 10⁻⁴. As an additional control against potential LBA artifacts that might result from the inclusion of DPANN, we also repeated the analysis, including inference of the underlying single-gene trees, without DPANN. The confidence set for this reduced analysis consisted of two trees, a tree placing the root between Euryarchaeota and the TACK/Lokiarchaeum lineage consistent with the full analysis, and a tree in which the root was placed on the branch leading to Lokiarchaeum (SI Appendix, Table S4). Although the Lokiarchaeum root could not be rejected in the reduced analysis, it was rejected by the full dataset, and it is not supported by our outgroup rooting analysis (SI Appendix, Fig. S13) or by published phylogenetic and comparative genomic analyses that group Lokiarchaeum with other “Asgard” Archaea within the TACK lineage (6, 7, 71). None of our analyses based on patterns of gene DTL provided any support for a root within a paraphyletic Euryarchaeota (10, 18).

DTL modeling approaches have been developed only recently, and their limitations are still being evaluated. To determine the robustness of our results, we performed a series of sensitivity and simulation analyses, which analyses indicated that our root inference is robust to high rates of horizontal transfer and variation in species sampling among gene families, and that our method robustly recovers the true root on simulated data (SI Appendix). An additional source of DTL error might be a kind of “small genome attraction,” in which the model favors a root that divides smaller from larger genomes on the tree. To investigate whether this might have been responsible for the support for basal DPANN, we repeated the rooting analysis using the 2,492 gene families that included at least one sequence from a DPANN archaeon. The 5% confidence set for the analysis of this reduced dataset contained only two rooted trees (SI Appendix, Table S5) and in both cases the root was placed between DPANN and all other Archaea. The difference between the two trees again lay in the position of the Thermoccocales, which were placed at the base of either the TACK Archaea or the Euryarchaeota. All of our sensitivity analyses agreed with the full data set in rejecting a root on Lokiarchaeum or among the Euryarchaeota.

Reconstructing Ancestral Archaeal Metabolisms.

Our DTL analysis provides an inference of the history of gene family evolution, including estimates of ancestral genome content (Fig. 3). To reconstruct ancestral metabolisms, we assigned functional annotations to the genes predicted to be present at each internal node on the tree, and mapped these onto core archaeal metabolic pathways (6). It is important to realize that these reconstructions are necessarily incomplete, because it is possible to reconstruct the history only of gene families that have survived to the present day in at least one of the sampled genomes. Moving back in time, the probability that genes on ancestral genomes survived to the present day decreases, and we estimate that 41% of the gene families that were present on the genome of the archaeal common ancestor have since gone extinct (Materials and Methods). These extinction probabilities can be used to correct ancestral genome size estimates (see below), although the functions of the extinct genes remain unknown.

Our analyses (Fig. 4) suggest that the LACA was an anaerobe that may have been able to fix CO₂ to acetyl-CoA via the Wood–Ljungdahl pathway and to subsequently generate acetate and ATP using an acetyl-CoA synthetase (arCOG01340) (SI Appendix). It is challenging to map genes to deep nodes in the phylogeny with high probability, and our reconstruction did not allow us to determine whether the electron donor for this reaction was organic or inorganic. Our mapping of the Wood–Ljungdahl pathway to the deepest nodes of the archaeal tree is in agreement with biochemical arguments and recent analyses using different methods (72, 73) that have suggested that the reduction of CO₂ with H₂ to produce organic compounds was central to the metabolism of an anaerobic last common ancestor of Bacteria and Archaea. Our analysis also suggests that the LACA had most of the modern archaeal transcription, translation, and DNA replication machineries, components of the exosome and proteasome, a secretion system, and some of the key genes for synthesizing archaeal ether lipids (SI Appendix, Figs. S16 and S17 and Table S6). The presence of lipids with a glycerol-1-phosphate backbone is often considered a hallmark of the Archaea and has been used to argue against the idea that eukaryotes originated from the endosymbiosis of the mitochondrial ancestor within an Archaeon (19, 20, 74, 75). Interestingly, we could not confidently map the origin of the key enzyme glycerol-1-phosphate dehydrogenase, which synthesizes the glycerol-1-phosphate backbone of canonical archaeal lipids, to the base of the Archaea, because it is absent from the published genomes or metagenomes of group II/III euryarchaeota, several members of the DPANN, and Lokiarchaeum (76). At least in the case of the Euryarchaeota and Lokiarchaeum, these absences likely reflect recent lineage-specific losses, because most Euryarchaeota, and all other members of the Asgard archaea, do encode glycerol-1-phosphate dehydrogenase (7).

Moving beyond the LACA, our analyses suggest that the Euryarchaeota/TACK common ancestor was also an anaerobe, possessing enzymes including superoxide reductase/desulfoferredoxin (pfam01880) commonly found in modern anaerobic and microaerophilic organisms. This ancestor also might have possessed an anaerobic proton-pumping system comprising membrane-bound F420- and/or H₂-dependent hydrogenases. Some of the recently discovered anaerobic Archaea that branch near the base of Euryarchaeota or TACK, such as Lokiarchaeum (6, 77), the Hadesarchaea (78), and some Bathyarchaeota (79), also have retained genes of the Wood–Ljungdahl pathway. Whereas some of the key enzymes of methanogenesis could be mapped (with P > 0.5) only to the base of the Euryarchaeota (Fig. 4), the recent discovery of methyl-CoM reductase in large-genome Bathyarchaeota (79) is consistent with an early origin of methane metabolism in Archaea (79), as is evidence of the presence of microbial methane—today produced exclusively by Archaea—in 3.46-billion-y-old rocks (19, 80).

Our analyses indicate that oxidative phosphorylation as attested by terminal oxidases and NADH dehydrogenase appears to have been acquired independently in several descendent lineages, including the TACK Archaea after their divergence from Lokiarchaeum and the stem leading to the Haloarchaeota (81). It is tempting to speculate that these parallel acquisitions of oxidative metabolisms may have been associated with the rise in atmospheric oxygen beginning around 2.5–2.3 billion y ago (82). Some of the genes today involved in sulfur metabolism also appeared first in the Euryarchaeota/TACK ancestor, including a potential sulfhydrogenase. Others, particularly genes for sulfur reduction, appear to have originated independently along the stems leading to different crenarchaeotal and euryarchaeotal lineages.

The inferred metabolic map of the DPANN common ancestor is similar to that of the LACA, consistent with the anoxic environments from which many members of this lineage have been obtained (5). However, in contrast to the LACA, the DPANN ancestor also encodes additional components of central metabolism, including enzymes involved in glucose and pentose sugar metabolism. The reconstruction suggests that the DPANN common ancestor may have been capable of anaerobic proton pumping via a V-type ATP synthase, given that two subunits of this membrane complex were mapped to the root of DPANN (SI Appendix, Table S6). Some modern DPANN species have subsequently lost these subunits, and it has been suggested (5, 83) that these may have a fermentative, parasitic, or symbiotic lifestyle.

Ancestral Growth Temperatures.

Previous work exploiting the correlation between sequence composition and optimal growth temperatures (OGTs) suggested that early Archaea were (hyper)thermophiles, with mesophily arising more recently in archaeal evolution (84, 85). Given that some DPANN genomes have been obtained from mesophilic environments, we investigated the impact of a basal DPANN clade on estimates of ancestral temperature. Our 45-gene alignment displayed a strong correlation between amino acid composition and OGT for modern Archaea (SI Appendix, Fig. S18), allowing us to infer temperature optima for ancestral nodes in the tree. We sampled 100 ancestral sequences for each node at the base of the tree using the branch-heterogeneous CoaLA model (84), performing the analysis both with and without DPANN (Materials and Methods). The LACA and the last common ancestors of each of the major archaeal clades (DPANN, Euryarchaeota+TACK/Lokiarchaeum, Euryarchaeota, and TACK+Lokiarchaeum) were all inferred to be thermophiles, and these inferences were robust to the inclusion of DPANN in the analysis (SI Appendix, Table S7); the median optimal growth temperature estimate for the LACA was 73.1 °C in the full analysis, and 75.7 °C in the analysis without DPANN. Interestingly, our model predicts mesophilic optimal growth temperatures for most modern DPANN genomes, consistent with the idea (84, 85)⁠ that adaptation to mesophily from a thermophilic ancestor occurred independently in each of the major archaeal clades.

Inferring Ancestral Genome Sizes.

The DTL model provides inferences of ancestral genome size, and, because the reconciliation model explicitly allows for horizontal transfer as well as gene loss, there is no trend toward inferring increasing genome size for earlier nodes on the tree. Thus, the use of this model ameliorates the “genome of Eden” (86) problem, a tendency toward inferring unrealistically large ancestral genomes in the absence of HGT that is so marked that it has been used to set a lower bound on rates of HGT through time (87). Previous simulation studies (67) and analyses of empirical data (88) have suggested that ancestral gene content inferences under this model are realistic and robust to gene tree uncertainty, and thus the ancestral sizes that we present here have been corrected to account for gene families that have been lost in all sampled species, as described above. Our analyses suggest that there has been an ongoing increase in gene content throughout archaeal history, from ∼1,090 genes in the common ancestor to 537–5,359 (mean, 1,686.4) genes among modern lineages. This trend is not dependent on the basal placement of the DPANN clade in the tree; in the analysis without DPANN, the common ancestor was predicted to encode 1,328 genes, increasing to 1,373–5,359 (mean, 2,081.1) genes among modern genomes. These reconstructions do not support the hypothesis, based on an analysis of phylogenetic presence-absence profiles (89), that a large-genome archaeal common ancestor gave rise to modern lineages by genomic streamlining.

Dynamics of Archaeal Genome Evolution: Gene Transfers, Duplications, and Losses.

Our reconciliations suggest that archaeal gene family evolution has been largely vertical (see also ref. 26), because for the majority (15,623) of gene families, vertical transmission events outnumber horizontal transfers [transfer ratio (TR) <0.5] (Materials and Methods). Interestingly, the distribution of TRs is multimodal, with a small peak of genes at TR >0.5 (SI Appendix, Fig. S19). In agreement with previous work on the transferability of genes with different kinds of functions (90), functional category had a significant effect on TR (P = 7.26 × 10⁻¹³⁴, ANOVA), with genes involved in carbohydrate metabolism (COG category G; P = 2.5 × 10⁻¹⁰, Fisher's exact test) and defense functions (COG category V) enriched in the set of frequently transferred genes with TR >0.5 (P = 6.7 × 10⁻¹², Fisher's exact test). Despite the overall predominance of vertical inheritance and the observed functional biases associated with HGTs, the cumulative effect of HGT on archaeal genomes is striking, and HGTs outnumber gene duplications on most (96 of 119) branches. Note that our inferences regarding HGT may represent underestimates, because increased taxon sampling may suggest that some inferred duplications instead are HGTs among close relatives. Remarkably, our reconstruction suggests that all of the gene families present at the root have experienced at least one HGT during archaeal evolution. Only 136 archaeal gene families are inferred to have entirely escaped HGT, and these are all recent originations.

In the discussion that follows, we define gene acquisitions as the sum of new genes that arise on a branch owing to lineage-specific innovation (i.e., apparently new genes with no detectable similarity to sequences from outside the subtended clade) or are obtained by HGT. The distributions of gene acquisition, duplication, and loss rates are continuous and correlated across the archaeal tree (all correlation coefficients ≥0.31; P < 0.01) (SI Appendix, Fig. S20), in agreement with some previous analyses of archaeal genome evolution (91). Acquisitions and duplications also are positively correlated with branch lengths (P < 0.05; SI Appendix, Fig. S20). Thus, according to our analyses, archaeal evolution is generally characterized by steady rather than punctate rates of genome change, with more events occurring on longer branches of the tree. Nonetheless, distributions for all of these processes have clear outliers (Fig. 5), indicating that some branches on the archaeal tree are exceptions to these background rates (SI Appendix, Table S8). The branch with the greatest number of gene acquisitions and duplications is that leading to the composite Lokiarchaeum genome. The high numbers may reflect the origin of some duplicated gene families shared with eukaryotes in the common ancestor of both lineages. However, the redundancy of the Lokiarchaeum composite genomic bin, which is estimated to include contigs from 1.4 closely related strains, also may be inflating estimates (6). The tip lineage with the second-highest number of acquisitions is that leading to Nitrososphaera gargensis, a group 1.1b Thaumarchaeote inhabiting a hot spring environment (92). This metabolically versatile archaeon has a much larger genome size (2.83 Mb) and gene complement (3,566 ORFs) than group 1.1a Thaumarchaeota, and has been inferred to have undergone extensive gene duplication, de novo gene origination, and horizontal acquisition (93). The crenarchaeote Sulfolobus solfataricus and the euryarchaeote Methanomassiliicoccus luminyensis also have experienced a large number of gene acquisitions, in both cases including a range of mineral and nutrient transporters.

The two stem lineages in which we observe the greatest number of gene acquisitions are the branches leading to the Haloarchaea and the Thaumarchaeota, two lineages that have undergone significant ecological transitions. Haloarchaea are suggested to have evolved into oxygen-respiring, light-harvesting heterotrophs from a methanogenic ancestor (81), whereas Thaumarchaeota may have evolved an ammonia-oxidizing lifestyle from an anaerobic ancestor (94). Horizontal transfer of metabolic genes from Bacteria has been implicated as an important process in these transitions (94–96), although the number of inferred transfers is sensitive to both the method used for mapping (26, 91, 94) and the taxonomic sampling of the lineages involved (97). Because the taxon sampling in our study was optimized for rooting the entire archaeal tree, our sampling within each of these groups is limited (two Haloarchaea, four Thaumarchaeota). Thus, from our analyses, it is difficult to pinpoint when these transfers occurred during haloarchaeal and thaumarchaeotal evolution.

Although we infer different numbers to published analyses, our findings are consistent with a substantial number of functionally relevant HGTs among Bacteria, Haloarchaea, and Thaumarchaeota. These include components of the electron transport chain and membrane transporters in Haloarchaea (SI Appendix, Table S9) and key components of the ammonia-oxidizing machinery in Thaumarchaeota (SI Appendix, Table S10), in agreement with recent large-scale phylogenetic analyses of genes shared between bacteria and mesophilic archaea (94–96). Interestingly, both stems also are inferred to have experienced a relatively large number of de novo gene originations and expansions of ancestral archaeal families. In the case of Haloarchaea, we identified 379 expanded or acquired gene families, including 109 (29%) corresponding to de novo gene origins, 156 (41%) corresponding to expansions of ancestral archaeal gene families, and 114 (30%) potential interdomain HGTs. These families have homologs in bacteria, although resolving the direction of transfer is difficult given the present data and methods. For Thaumarchaeota, we identified 17 de novo origins (16%), 72 expansions (69%), and 15 interdomain HGTs (15%). The haloarchaeal stem was the branch experiencing the greatest number of gene losses.

Conclusion.

In the present study, we used large amounts of genomic data and a method that implicitly considers patterns of genome duplication, HGT, and gene loss (67) to generate a rooted tree for the Archaea, one of the two primary domains of life (19). The DTL model performed well in simulations and in our case used phylogenetic signals from 31,236 homologous gene families, compared with the small universal core of single-copy orthologous genes typically used for outgroup rooting. The DTL analyses infer a new root between the DPANN clade and all other Archaea, with the Euryarchaeota and the TACK/Lokiarchaeum clade resolved as monophyletic sister lineages. Monophyly of DPANN was supported by supertrees, supermatrices, and DTL modeling, and thus, notwithstanding legitimate concerns about potential LBA artifacts, is the hypothesis best supported by our analyses. Its robustness will be tested as methods and genomic sampling of the relevant groups continue to improve.

The DTL analysis and new root provides inferences of gene content evolution that are consistent with inferences of early archaeal physiology based on other lines of evidence (20, 72, 73). Our analysis suggests that the LACA was an anaerobe that fixed carbon via the Wood–Ljungdahl pathway, and that adaptations to aerobic metabolism evolved independently across the tree. We infer that ecological transitions within the Archaea are associated with substantial gene content turnover, involving both HGT and the evolution of lineage-specific genes. Although our analyses agree that HGT is an important feature of archaeal evolution, the majority of transmission events appear to be vertical rather than horizontal, preserving a strong vertical trace between lineages. In contrast to hypotheses in which a large-genomed archaeal common ancestor gave rise to modern lineages by streamlining (89), the DTL analyses imply a moderate increase in gene content throughout archaeal history from a common ancestor that had a relatively small genome. Our analyses also suggest that adaptation to mesophily from a thermophilic ancestor occurred independently in each of the major archaeal clades.

Sequences and Alignments.

We used the OMA algorithm (22) to identify orthologous gene families on 62 archaeal genomes, resulting in 4,664 orthologous families with at least four members. Families were screened for interdomain HGT using a BLASTP-based protocol requiring that all genes be more similar (lower E-value) to one or more sequences from other archaeal genomes than to any sequences from bacteria or eukaryotes. Sequences that did not meet this requirement were filtered out, resulting in a set of 3,266 orthologous archaeal protein families containing four or more sequences. Sequences were aligned using MUSCLE (100), and poorly aligning regions were identified and removed using BMGE (47) under the BLOSUM30 substitution matrix.

Single Gene Trees.

Single gene trees were inferred using the C60+LG model in PhyloBayes 4.1 (101); this model is optimized for smaller datasets on which more general models (102) can show convergence problems. Two chains were run in parallel, and convergence was assessed using the bpcomp and tracecomp programs in PhyloBayes. A consensus tree was built once the maximum interchain discrepencies in bipartition frequencies and a range of continuous model parameters had dropped below 0.1, with effective sample sizes for continuous parameters >100. One-quarter of sampled points were discarded as burn-in.

Supermatrices.

We identified 57 orthologous families that were present on all, or all but one, of the completely sequenced genomes in our dataset. This slightly relaxed criterion allowed for occasional gene losses or misannotations on published genomes. Because assemblies for uncultured organisms are often incomplete, we did not require the presence of these marker genes on all single-cell genomes, but did include orthologs where present. The gene trees from these families were visually inspected, and only those recovering the monophyly of both the TACK Archaea and the core Euryarchaeota at PP >0.7 were concatenated. This resulted in a final set of 45 single-copy orthologous marker genes and a concatenation of 10,738 aligned amino acid positions. The LG, CAT+Poisson, and CAT+GTR substitution models were fit to this concatenation using PhyloBayes-MPI (103), with posterior predictive simulations used to evaluate model adequacy. None of these models was able to adequately account for the across-site (z = 6.14; P = 0, CAT+GTR) or across-branch compositional heterogeneity (z = 8.001; P = 0, CAT+GTR), potentially leading to phylogenetic artifacts. Thus, we explored data-recoding techniques as a means of ameliorating these compositional biases. Even after data recoding, the data contained both across-site and across-branch compositional heterogeneity that was not adequately anticipated by the model (P = 0, CAT+GTR+Dayhoff4), but a reduction in the z-scores associated with the posterior predictive tests (z = 3.43 for across-site compositional heterogeneity, z = 4.66 for across-branch compositional heterogeneity) suggested improved model fit.

Supertrees.

We used MRP (25) to infer a supertree from majority-rule posterior consensus trees for the orthologous archaeal protein families. Out of 3,266 families, 24 produced comb trees and 3,242 had at least some resolution at PP >0.5; we used this latter set in the MRP analysis. The input trees were definitive, producing a single most parsimonious supertree.

Modeling Gene DTL.

We built an expanded set of gene family trees that included both paralogs and orthologs. We performed an all-versus-all BLASTP analysis of our 62 archaeal genomes (E-value threshold <10⁻⁵), then inferred gene clusters using MCL (Markov clustering algorithm) with an inflation parameter of 1.4. Sequence sets were aligned and masked, and gene trees were inferred as above. Because these homologous gene families were larger than the orthologous gene families used above, we used slightly less stringent convergence criteria to obtain gene tree samples within a tractable amount of time (maximum difference in bipartition frequencies, 0.3; minimum effective sample sizes of continuous parameters, >50). We performed gene tree-species tree reconciliation using the ALEml_undated algorithm of the ALE package (67) (https://github.com/ssolo/ALE), which uses a probabilistic approach to exhaustively explore all reconciled gene trees that can be amalgamated as a combination of clades observed in a sample of gene trees. We estimated a single global set of ML DTL rates for each rooted species tree. In estimates of ancestral gene content, we used extinction probabilities conditional on the estimated ML rates to account for genes that have gone extinct.

Functional Annotation of Gene Families.

ArCOG categories were assigned by BLASTing each member of each gene family against the 2015 version of the ArCOG database (104). For KEGG Orthology number assignment, we selected the medoid (the sequence with the shortest summed genetic distances to all other sequences in the family, calculated under the BLOSUM62 substitution matrix) for annotation, and searched this against the KEGG database (105) using the KAAS annotation server (106) (accessed February 2016). We reconstructed ancestral gene family repertoires from the DTL model by selecting all families predicted to be present at a given node with P ≥ 0.5. We assessed the metabolic capabilities of ancestral genomes using the KEGG Module tool. Genes gained along specific branches of the archaeal tree were identified by screening for gene families whose size increased by ≥1 along that branch. The origins of these genes were assessed by BLASTing against our sample of archaeal genomes and a reference set of bacterial genomes. Genes gained on a branch with significant hits to one or more bacterial genomes (BLASTP E-value <0.00001) but no homology to any other archaeal genomes were classified as putative bacterial HGTs, genes with no homology to any other sequenced genome were classified as de novo gene originations, and genes with homology to other archaeal genomes were classified as expansions of ancestral archaeal families, whether by gene duplication or by within-Archaea HGT transfer.

Ancestral Temperature Estimates.

Estimates of ancestral OGTs are based on reconstructed ancestral sequences based on three different concatenates, analyzed under the LG+4G+COaLA model. At the root, 19 free parameters were estimated for the amino acid frequencies. Parameters were estimated by ML using the Bio++ libraries (107). A correspondence analysis on amino acid composition was performed for each of the three concatenates using the R package ADE4 (108). In all cases, axis 2 was found to correlate with OGT. The OGT for ancestral nodes was predicted by linear regression (predict.lm in R, with “interval= confidence”) to take uncertainty in the parameter estimates of the regression line into account.