Fieldwork, sampling, and DNA extraction of lava rock samples from active fumaroles in Mauna Ulu volcanic area in Hawai‘i were previously reported, and the MAGs in this study originated from sample numbered 86B (Cockell et al., 2019; Hughes et al., 2019). Biofilm samples from two hot springs in El Tatio geothermal field in Chile were collected during a field trip in August 2018. The MAGs reported here came from samples collected at two different hot spring features: Cacao East 5 (abbreviated as CE5; GPS coordinates of −22.350367, −68.008050) and Poppy (GPS coordinates of −22.333478, −68.013011) pools. A field photo of Poppy pool is shown in Supplementary Figure S1 and that of CE5 has been reported previously (Megevand et al., 2022). At the time of collection, the main CE5 pool had a temperature of 83.4°C and pH of 7.2, whereas the main Poppy pool had a temperature of 78.3°C and pH of 6.85. Sample material for DNA extraction from Poppy consisted of a brown-colored biofilm scraped from a piece of hot spring sinter that was fractured with sterile tweezers off the rim of the actively splashing hot spring. Sample material for DNA extraction from the more quiescent CE5 pool was collected with a sterile tweezer as a vertical section of a completely submerged hot spring mat located just below the water line at the pool’s edge. After collection, the hot spring sinter and mat samples were placed into sterile Qorpak jars (Berlin Packaging, Clinton, PA), and transported within a few hours to a hotel where they were stored at ~4°C overnight. The sample jars were then transported on ice to the PNNL laboratory where they were stored at −80°C until they were shipped on blue ice to the University of Hawai‘i laboratory, where they remained frozen until they were processed for DNA extraction. DNA was extracted from 0.5 g of each sample using the Qiagen DNeasy PowerSoil DNA kit and quantified with a Qubit instrument.

Illumina libraries of samples from Chile were prepared using NEBNext Ultra DNA II Library Preparation Kit (New England Biolabs, Cat. #E7645L). Sequencing libraries were generated from sample material that had input DNA amount varying between 5 and 100 ng. The DNA was fragmented with a Covaris E220, the ends made blunt, and adapters and indexes added to blunt ends for sequencing on an Illumina sequencer. Illumina libraries were eluted in DNA Elution Buffer (Zymo Research, Cat. #D3004-4-10), and concentrations of the libraries were obtained using the Qubit dsDNA HS Assay (ThermoFisher Scientific, Cat. #Q32854). The average size of the DNA in the library was determined by the Agilent High Sensitivity DNA Kit (Agilent, Cat. #5067–4,626). Libraries were quantified with the Library Quantification Kit – Illumina/Universal Kit (KAPA Biosystems, KK4824), and sequenced using a NextSeq v2.5 Reagent Cartridge (Illumina, Cat. #20024908).

Raw metagenomic reads were first trimmed and filtered for contamination using BBTools v38.87 (Bushnell, 2014). Quality filtering excluded adapter contaminants, read regions with scores < 20, and all reads under 50 bp. The resulting paired metagenomic reads were assembled using MetaSPAdes v3.14.0 (Nurk et al., 2017) with kmers specified as 21, 33, 55, and 77. MetaBAT2 v2.15 (Kang et al., 2019) was used to create genome bins from the assembled contigs by first using Seqtk v1.3 (Li, 2018) to remove contigs under 1 kb, then mapping the filtered contigs using BBTools v38.87 (Bushnell, 2014) with default parameters. Samtools v1.10 (Li et al., 2009) was used to sort and index the mapped contigs using default parameters, and BBTools v38.87 (Bushnell, 2014) then generated contig depths for all contigs over 1.5 kb. This contig depth file containing mean and variance of the base coverage depth was used in generating genome bins with MetaBAT2 v2.15 (Kang et al., 2019) with otherwise default parameters. The lineage workflow parameter of CheckM v1.1.3 (Parks et al., 2015) was used to assess completeness and contamination of the binning results to ensure bin quality for downstream analysis.

GTDB-Tk v1.4.0 (Chaumeil et al., 2020) was used to assign taxonomic affiliation to all genome binning results. All bins identified by GTDB-Tk as Caldarchaeales were selected for downstream analysis. QUAST v5.0.2 (Gurevich et al., 2013) was used to determine the number of scaffolds, genome size, N50 value, and genomic G + C content. Estimated genome size was determined from CheckM (Parks et al., 2015) completeness and the genome size result from the QUAST analysis. FastANI (Jain et al., 2018) was used to compute ANI values between the five Caldarchaeales MAGs and 79 Caldarchaeales genomes downloaded from the GTDB release 207_v2.

Prodigal (Hyatt et al., 2010), with the “-p meta” option, was used to predict protein-coding genes in the five Caldarchaeales MAGs. All predicted proteins were analyzed with InterProScan (Jones et al., 2014) with default parameters to annotate protein domains and assigned to archaeal clusters of orthologous genes (arCOGs; Makarova et al., 2015b) by eggNOG-mapper v2.1.9 (Cantalapiedra et al., 2021) with default settings. Protein-coding genes were also annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000) tool GhostKOALA v2.2 (Kanehisa et al., 2016) and METABOLIC (Zhou et al., 2022); and annotations for genes of interest were confirmed by the BLASTp (Altschul et al., 1997, 2005) searches that were performed against the following databases: NCBI-nr (Sayers et al., 2021), UniRef100 (Suzek et al., 2015), and UniProtKB reference proteomes plus Swiss-Prot (UniProt Consortium, 2021; E-value cutoff ≤ 1e⁻⁵). Subcellular localization of peptidases found by METABOLIC against the MEROPS database was predicted by the PSORTb web tool v3.0.3 (Yu et al., 2010). The number of rRNA-coding sequences and presence of 16S rRNA gene in the MAGs were determined with Anvi’o (Eren et al., 2021). The Hidden Markov Model (HMM) profiling function, which utilizes Barrnap,¹ predicted the locations of rRNA genes. The tRNAscan-SE software v1.3.1 (Chan and Lowe, 2019) was used to determine the number of tRNAs in each MAG. Metabolic pathways were reconstructed based on the hits for key metabolic marker genes and the reference pathways depicted in KEGG (Kanehisa and Goto, 2000) and MetaCyc (Caspi et al., 2020) databases.

CRISPR arrays of repeat-spacer units, and cas genes were identified using the combination of web tools CRISPRCasFinder (Couvin et al., 2018), CRISPRone (Zhang and Ye, 2017), CRISPRminer2 (Zhang et al., 2018), and CRISPRCasTyper (Russel et al., 2020). CRISPRone, CRISPRminer2, and CRISPRCasTyper were run with default settings. While running CRISPRCasFinder, the “General” clustering model and the “Unordered” function were used to identify cas genes. The CRISPRCasFinder tool has an evidence level rating that ranges from 1 to 4 for detection of CRISPR arrays. Only CRISPR arrays having evidence levels 3 and 4 were included in our analysis. While defining the cas locus, CRISPR locus, and CRISPR-Cas locus, and determining the type of a cas locus here, we followed published criteria (Zhang and Ye, 2017): (1) a cas locus should contain at least three cas genes, and at least one of those should belong to the universal cas genes for CRISPR adaptation (cas1 and cas2) or the main components of interference module including cas7, cas5, cas8, cas10, csf1, cas9, cpf1 (Makarova et al., 2015a); (2) CRISPR arrays that are close to each other (≤200 bps) and share very similar repeat sequences were considered to be in the same locus; and (3) the CRISPR(s), together with nearby (within 10,000 bps) cas genes, are defined as a CRISPR-Cas locus. The type of each cas locus was determined based on type signature cas genes listed in Makarova et al., 2015a,b The type-assignment of a cas locus was considered confident if it had at least three type-consistent signature cas genes, except for type V, which has only one signature gene, cpf1 (Makarova et al., 2015a; Zhang and Ye, 2017). If a cas locus included only one or two type signature cas genes, we considered the type-assignment as putative.

The VIBRANT v1.2.0 (Kieft et al., 2020) with default settings through the CyVerse Discovery Environment at https://de.cyverse.org/de (Merchant et al., 2016) and VirSorter v1.0.5 (Roux et al., 2015) with default settings against the Viromes database through the DOE Systems Biology KBase at http://kbase.us (Arkin et al., 2018), were used on all contigs of at least 1 kb to identify putative viral genes in the MAGs. All putative viral gene sequences detected by VIBRANT and VirSorter were queried against UniRef100 (Suzek et al., 2015), and UniProtKB reference proteomes plus Swiss-Prot (UniProt Consortium, 2021) databases using the BLAST tool (Altschul et al., 1990) with BLASTp (Altschul et al., 1997, 2005) at https://beta.uniprot.org/blast for functional annotation. Furthermore, these putative viral gene sequences were also functionally annotated by querying them against the NCBI non-redundant protein sequences (nr; Sayers et al., 2021) and IMG/VR Viral Resources v3 (Roux et al., 2021) databases using BLASTp. The BLAST hits with e-values equal to or lower than 1e⁻⁵ were accepted as statistically significant hits. HHpred, HHblits, and ProtBLAST/PSI-BLAST tools of the MPI Bioinformatics Toolkit (Zimmermann et al., 2018; Gabler et al., 2020) and NCBI CD (Conserved Domains)-Search tool (Lu et al., 2020) were also utilized for further annotation of putative viral genes in our MAGs. The search for auxiliary metabolic genes (AMGs) in the putative prophage regions of the MAGs was carried out based on the list of AMGs provided by Kieft et al. (2020).

VIBRANT uses HMM profiles from three databases: KEGG KoFam, March 2019 release (Kanehisa and Goto, 2000; Aramaki et al., 2020), Pfam (v32; Mistry et al., 2021), and Virus Orthologous Groups (VOG) release 94² to identify and annotate viral proteins (Kieft et al., 2020). On the other hand, VirSorter compares predicted protein sequences to Pfam (v27; Mistry et al., 2021) and to the user-selected viral reference database (Viromes database was chosen for this study) to detect viral-like genes (Roux et al., 2015).

If the protein sequence of a gene found in a putative viral region did not have any hits in KEGG, KoFam, Pfam, VOG, and Viromes databases, and this protein sequence only yielded hits to hypothetical or uncharacterized proteins when BLASTed against UniRef100, UniProtKB reference proteomes plus Swiss-Prot, NCBI nr, and IMG/VR databases, it was then annotated as a hypothetical/uncharacterized protein. We also manually inspected the KEGG, eggNOG-mapper, and InterPro annotations of the MAGs to search for viral hallmark proteins, specifically to identify any viral region not detected by VIBRANT and VirSorter tools. Viral hallmark proteins are necessary for productive infection, such as structural (e.g., capsid, tail, baseplate, portal, coat, spike), terminase, or viral holin/lysin proteins, and their presence is one of the metrics used to capture viral signatures in metagenomes (Roux et al., 2015; Kieft et al., 2020).

Pangenomic analyses were performed using Anvi’o v7.1 for three of the MAGs with sufficient genome completeness (>55%): 146_bin.21, 146_bin.25, and 1054_113_bin.10 and their closest relatives in the GTDB release 207_v2 (Table 1). Since two of the MAGs generated from the Chile samples had genome completeness levels lower than 55%, they were excluded from the pangenomic analyses. All of the closest relatives of 146_bin.21, 146_bin.25, and 1054_113_bin.10 in the GTDB displaying completeness levels higher than 55% were included in the analyses. The first pangenomic analysis included 1054_113_bin.10 and three reference genomes that showed >95% ANI (Average Nucleotide Identity) with 1054_113_bin.10 and belonged to the GTDB species Caldarchaeum subterraneum_B (GCA_011373365.1, GCA_011364015.1, GCA_011364605.1).

The second pangenomic analysis included 146_bin.25 and two reference genomes that showed >95% ANI with 146_bin.25 and belonged to the GTDB species JGI-OTU-1 sp011364265 (GCA_011364265.1, GCA_011369755.1). Lastly, the third pangenomic analysis included 146_bin.21 and two reference genomes that showed approximately 80% ANI with 146_bin.21 and belonged to the GTDB species WAQC01 sp015520425 (GCA_015520425.1, GCA_015522085.1). Due to relative incompleteness of closely related MAGs in GTDB, we were only able to use very few genomes in the pangenomic comparisons. To ensure continuity, previous protein-coding gene annotations completed by Prodigal were first processed using “anvi-script-process-genbank” and imported into the Anvi’o contigs databases using the “—external-gene-calls” flag and “anvi-import-functions” program. After the pan-genomes were generated using the “anvi-pan-genome” program, the ANIb values between query and reference genomes were imported using PyANI from the “anvi-compute-genome-similarity” program within the Anvi’o suite.

The Genome Taxonomy Database Toolkit (GTDB-Tk) version 2.1.1 was used alongside database release 207_v2 to classify the five MAGs presented here. The resulting multiple sequence alignment (MSA) containing 53 concatenated and conserved single-copy archaeal marker genes (10,153 characters) in 3492 archaeal genomes, which included the five MAGs generated in this study, was used as an input to perform a model test in the IQ-Tree program (version 2.1.4-beta; Minh et al., 2020). A maximum likelihood phylogenomic tree was constructed using the best model chosen based on the AIC and BIC criteria (LG + F + G4), with 1,000 ultra-rapid bootstrap replicates. The resulting tree was uploaded to iTOL (Interactive Tree of Life; Letunic and Bork, 2021) to collapse known clades, and edited in Inkscape for esthetics only.

Barrnap v0.9³ was used to identify and extract 16S rRNA gene sequences from three of our five Caldarchaeales genomes (1054_113_bin.10, 146_bin.21, 146_bin.25). To detect the presence of Caldarchaeales across diverse habitats, all three 16S rRNA gene sequences were used as queries for the Integrated Microbial NGS Platform (IMNGS; Lagkouvardos et al., 2016) to perform a Paraller similarity search (97% similarity), which analyzes all available prokaryotic 16S rRNA amplicon experimental datasets from the International Nucleotide Sequence Database Collaboration (INSDC). The INSDC includes datasets submitted from the DNA Data Bank of Japan (DDBJ), European Molecular Biology Laboratory’s European Bioinformatics Institute (EMBL-EBI), and National Center for Biotechnology Information (NCBI). The IMNGS result provided the relative abundance of reads and environmental descriptors for each hit, but data for location coordinates were searched manually in the NCBI SRA database or original literature due to inconsistencies in location reporting. The relative abundance was calculated by dividing the number of sequences identified by IMNGS to be >97% identical to the query 16S rRNA gene sequences by the total number of sequences in each 16S rRNA amplicon survey. All non-environmental 16S rRNA amplicon surveys (e.g., oral, skin, gut) were removed from downstream analysis, and amplicon sequences with identical location coordinate data were combined. The 16S rRNA gene sequences of our three Caldarchaeales MAGs were also searched against the Silva 16S rRNA database to obtain GenBank accession numbers for all hits. These GenBank accession numbers were downloaded and parsed using custom Python scripts to extract environmental information and coordinate location data. GenBank accessions containing location coordinate data were used in the geographic analysis, and again, identical locations were combined. Environmental features and location coordinates data for each hit were visualized in R using the ggplot2 (Wickham, 2016) and maps (Becker et al., 2022) packages. Relative abundance data were visualized for all 16S rRNA amplicon surveys found using IMNGS.

Searches against the IMNGS and Silva 16S rRNA databases revealed that the Caldarchaeales are represented in a diverse array of features, including hot springs, marine geothermal features, lake and ocean sediments, soil, and other aquatic systems. They were also determined to be present in low relative abundance (<0.0001%) in plants such as Populus deltoides, fish gut metagenomic studies, and in surveys of human skin, oral, and gut microbial communities. Although they appear to inhabit a wide array of ecological niches, Caldarchaeales are especially prevalent in hot springs and other geothermal features (Figure 1; Supplementary Table S1). We cannot exclude the possibility that some of the unusual habitats the Caldarchaeales are found might be attributed to barcode misassignments or contaminations in the samples we obtained from these public databases. The Great Basin and Yellowstone regions of the United States host the greatest relative abundance of these archaea. These regions likely provide optimal conditions for thermophilic archaea since the abundance of affiliated sequences in all other features and areas is at least 100-fold less than in US hot springs. It should also be noted that in some soil and sediment environments, temperatures can be similar to those in many hot springs, such as a fumarolic soil sample from Antarctica in which the temperature was 65°C (Herbold et al., 2014). However, there is a possibility that primer bias could have also contributed to higher abundances of these archaea in some of the samples. Additional sampling and metagenomic studies will provide more insights into the geographical distribution of the Caldarchaeales.

General features of the five MAGs along with genome completeness and contamination estimates are shown in Table 1. GTDB-Tk classification of the MAGs resulted in one from Hawai‘i (146_bin.21) being assigned species-level novelty in the GTDB-designated genus WAQC01, in the family Caldarchaeaceae, order Caldarchaeales. In the GTDB release 207_v2, only two genomes represent WAQC01, and these are MAGs generated from hydrothermal deposits in the subsurface of a submarine volcano in the Pacific Ocean (Reysenbach et al., 2020). The other MAG from Hawai‘i (146_bin.25) was assigned to the GTDB species JGI-OTU-1 sp011364265 in the GTDB genus JGI-OTU-1, family HR02, order Caldarchaeales. It forms a sister-level relationship with Candidatus “Calditenuis aerorheumensis.” Three MAGs from Chile were assigned to the GTDB species Caldarchaeum subterraneum B in the family Caldarchaeaceae, order Caldarchaeales.

IQ-Tree phylogenomic inference of the five MAGs alongside the archaeal backbone tree from the GTDB release 207_v2, shows the three bins from Chile in a clade with Ca. C. subterraneum (Figure 2). 146.bin.25 groups closely with members of family HR02 and forms a sister-level relationship with a SAG from Great Boiling Spring (Rinke et al., 2013). It is interesting to note that 146_bin.21 forms a monophyletic clade with two deep-sea hydrothermal vent-associated lineages (Reysenbach et al., 2020) despite it coming from a terrestrial environment. More importantly, we obtained a monophyletic clade comprising all Caldarchaeales as sister to the clade comprising the orders Geothermarchaeales (previously Candidatus Geothermarchaeota), Conexivisphaerales and Nitrososphaerales (formerly classified as the phylum Thaumarchaeota).

Based on the SeqCode guidelines to assign valid names to microbial genomes with relatively high completeness and low contamination, both MAGs from Hawai‘i can be given valid genus and species names as both are above 90% completeness, below 5% contamination, and have small numbers of contigs (Hedlund et al., 2022; Whitman et al., 2022). Here, we propose the names Pelearchaeum maunauluense^Ts gen. Nov., sp. nov. for 146_bin.21 after Pele, the goddess of volcanoes and fire, and Calditenuis fumarioli ^Ts sp. nov. for 146_bin.25, which belongs to the same genus as Candidatus “Calditenuis aerorheumensis” (ANI value of 81.07% between them; Figure 2). Due to monophyletic grouping of 146_bin.25 with 10 other members (Figure 2), we also propose to name the clade (GTDB family designation HR02) as Calditenuaceae and designate Calditenuis fumarioli^Ts as type species. A list of proposed names, etymologies, and nomenclatural types are shown in Table 2. Hereafter, we will refer to 146_bin.21 as P. maunauluense and 146_bin.25 as C. fumarioli.

Metabolic features inferred from the genomes of the five MAGs are summarized in Figure 3; Supplementary Tables S2, S3, S10. GhostKOALA and eggNOG-mapper annotations of all coding sequences can be found in Supplementary Tables S4, S5, respectively. Their metabolic potentials are discussed in greater details below.

The presence of some key marker genes for the 3HP/4HB and DC/4HB cycles in the MAGs from Hawai‘i hints at the capacity to fix carbon. However, these cycles are not complete. All MAGs from Hawai‘i and Chile also possess metabolic potential for roTCA cycle, but only further experimental validation will determine their functionality. Metabolic pathway predictions suggest the potential denitrification and tetrathionate reduction pathways might enable these lineages to respire diverse compounds other than oxygen. All MAGs except 1054_113_bin.10 encode reverse gyrase, a marker gene for hyperthermophiles, which is consistent with the environmental conditions of the habitats from which the MAGs were recovered. The occurrence of genes involved in aerobic and anaerobic respiration in the MAGs indicates a facultative anaerobic lifestyle, which would enable growth in anoxic environments. We also determined that the MAGs contain several genes implicated in oxidative stress tolerance, such as those for superoxide dismutase, thioredoxin reductase, and peroxiredoxin, which could be used to alleviate oxidative damage.

Some of our MAGs also possess genes involved in arsenic detoxification and resistance to mercury (Figure 3). Only P. maunauluense genome encodes genes capable of encoding arsenate reductase (arsC), though all other MAGs encode the gene for the ArsR family transcriptional regulator (arsR). Evidence for genes that encode the arsenical pump membrane protein (arsB) and arsenite/tail-anchored protein-transporting ATPase (arsA) were detected in all MAGs except 1054_113_bin.10. As for resistance to mercury, P. maunauluense and C. fumarioli genomes encode the gene for mercuric reductase (merA). The presence of arsenic and mercury in fumarolic rocks of Hawai‘i is currently unknown. Hot springs are known to contain various levels of arsenic and at El Tatio, the concentrations can sometimes reach 40–50 mg/L of hot spring water, which is among one of the highest concentrations measured in any natural habitat (Wang et al., 2018) and therefore presence of transporters to pump out various arsenic related compounds such as arsenous acid [As(OH)₃] appear to indicate that Caldarchaeales in these habitats may be actively utilizing these detoxification systems for their survival.

A recent study has reported the enrichment and stable laboratory growth of Caldarchaeales from sediments of Great Boiling Spring (Nevada, United States), designated Wolframiiraptor gerlachensis, which belongs to the GTDB family NZ13-MGT (Buessecker et al., 2022). Cultivation experiments and the analysis of a MAG from W. gerlachensis, which includes six annotated tungsten (W)-dependent ferredoxin oxidoreductases that could play central roles in anaerobic carbohydrate degradation, have suggested that tungsten could be essential for carbohydrate metabolism, and that it may play a critical role in one or more of the annotated W-dependent ferredoxin oxidoreductases, thereby indicating that the growth of W. gerlachensis requires the presence of tungsten, a biologically rare trace metal (Buessecker et al., 2022). In addition, the evolutionary analysis of 78 MAGs representing four genera and 11 species of the GTDB family NZ13-MGT, designated as Wolframiiraptoraceae, has demonstrated that putative tungstate transporters and tungstoenzymes were ancestral in this lineage. Furthermore, homologs of 3 W-dependent ferredoxin oxidoreductases: aldehyde:ferredoxin oxidoreductases (AOR), glyceraldehyde-3-phosphate:ferredoxin oxidoreductases (GAPOR), and formaldehyde:ferredoxin oxidoreductases (FOR) are widely distributed in the family Wolframiiraptoraceae genomes, and sparse in other Caldarchaeales.

Also, the majority of the Wolframiiraptoraceae genomes has been detected to encode Tup ABC tungstate transporter system, whereas other ABC tungstate transporter systems: Wtp and Mod have only been identified in two species in the family Wolframiiraptoraceae (Buessecker et al., 2022). Given the significant role of tungsten in sugar and peptide catabolism within the Caldarchaeales (Buessecker et al., 2022), we searched for tungstate transporters and tungsten-associated enzymes in the five MAGs. Five AOR homologs were identified, one in 1054_113_bin.10, and four genes in P. maunauluense, three of which are fragmented. However, we did not recover any GAPOR or FOR homologs in the MAGs. Pelearchaeum maunauluense genome also encodes three subunits of the Wtp molybdate/tungstate transporter system, but no tungstate transporter system was annotated in 1054_113_bin.10, possibly due to incomplete genomes (Figure 3). Intriguingly, C. fumarioli includes three subunits of the Wtp molybdate/tungstate transporter system and an adjacent gene encoding a molybdate transport system regulatory protein, but no W-dependent ferredoxin oxidoreductases. The AOR homologs encoded in P. maunauluense and 1054_113_bin.10 may be involved in the protein catabolism by converting aldehydes, which were produced from 2-Keto acids by ferredoxin oxidoreductases to organic acids and also in the Branched Entner-Doudoroff pathway; GAPOR is only known to function in unusual Embden-Meyerhof pathways in which it converts glyceraldehyde-3-phosphate to 3-phosphoglycerate, replacing glyceraldehyde-3 phosphate dehydrogenase and phosphoglycerate kinase enzymes used in the canonical Embden-Meyerhof pathway (Mukund and Adams, 1995; Bevers et al., 2005). Pelearchaeum maunauluense, C. fumarioli, 1054_108_bin.3, and 1053_113_bin.10 encode both glyceraldehyde-3 phosphate dehydrogenase and phosphoglycerate kinase enzymes, which might explain the absence of a GAPOR homolog from these genomes.

Previously, integrated mobile genetic elements (IMGEs) that include bacterial insertion sequence (IS)-like transposons, prophages, and cryptic integrated elements were reported in several Caldarchaeales genomes (Nunoura et al., 2011; Hua et al., 2018; Krupovic et al., 2019). In line with these studies, we also found putative prophage regions and genes encoding transposases in our MAGs, suggesting the MAGs might harbor IMGEs. Prophage-related genes and transposases identified in the MAGs are shown in Supplementary Tables S6, S7, respectively.

Viral genes were not detected in the MAG 1054_113_bin.10. However, the MAGs of P. maunauluense, C. fumarioli and 1054_113_bin.2 include prophage regions, and MAG 1054_108_bin.3 has three viral contigs, which we refer to as putative prophage regions (Figure 4). Although VIBRANT and VirSorter predicted the three contigs of MAG 1054_108_bin.3 to be entirely viral, they very likely belong to prophage regions integrated into the host genome, given MAG 1054_108_bin.3 is a very fragmented and incomplete genome of 167 contigs and 44.34% completeness.

Prophage regions in the MAGs comprise numerous viral hallmark proteins such as the DNA packaging enzyme terminase, phage integrase, a prohead protease involved in phage assembly and maturation, and structural proteins including phage capsid, tail, and portal proteins. In addition to viral hallmark proteins, prophage regions encode proteins associated with viral genome replication such as Staphylococcal nuclease homolog, which is thermostable, replicative DNA helicase Mcm (minichromosome maintenance helicase), RNA 2′,3′-cyclic 3′-phosphodiesterase, and ParB family transcriptional regulator. The replicative DNA helicase Mcm was shown to be frequently recruited from the host as the main replication protein of various crenarchaeal and euryarchaeal mobile genetic elements, including viruses and plasmids (Krupovic et al., 2010, 2019; Kazlauskas et al., 2016). RNA 2′,3′-cyclic 3′-phosphodiesterase was also shown to function as a highly heat-stable cyclic nucleotide phosphodiesterase with GTP-dependent RNA ligase activity in the hyperthermophilic archaeon Pyrococcus furiosus (Kanai et al., 2009). The presence of a thermostable nuclease and RNA 2′,3′-cyclic 3′-phosphodiesterase in the provirus regions suggests viruses infecting Caldarchaeales might have evolved strategies to replicate their genomes in high-temperature environments. The prophage regions annotated here also carry genes encoding functionally diverse proteins that might contribute to the growth of these archaea in extreme conditions. They are listed in Table 3. These findings are consistent with a previous study reporting that IMGEs found in 20 Thaumarchaeota genomes and one Caldarchaeales genome encode multiple AMGs, various dehydrogenases, stress response proteins, membrane transporters of cations and drugs, and chemotaxis protein receiver domains, which could contribute to the fitness, adaptation, and survival of the archaeal hosts by providing mechanisms that respond to stressful environmental conditions, and improve the host metabolic potential (Krupovic et al., 2019).

Viral gene clusters found in the MAGs might belong to Mu-like prophages, as the terminase enzymes encoded in the provirus regions of P. maunauluense, 1054_113_bin.2, and 1054_108_bin.3 share significant amino acid identities with Mu phage terminases, and one of the phage tail proteins encoded in the provirus region of C. fumarioli shares significant amino acid identities with that of the prophage Mu tail protein. The finding that some viral hallmark proteins in the provirus regions of the MAGs display homology to Mu phage proteins suggests that Mu-like phages could be infecting Caldarchaeales hosts. This finding also corroborates a previous study in which Mu-like prophages were found in two Caldarchaeales MAGs (Hua et al., 2018). Viruses identified in these MAGs might enhance the metabolic potential of their hosts through horizontal gene transfer events, facilitate niche-specialization, and may reduce competition for resources in extreme habitats such as hot springs and fumaroles found in the study sites presented here. Archaeal viruses are among the most enigmatic viruses known due to a relatively small number of them described to date compared to bacteriophages and eukaryotic viruses and due to their extraordinary genomic diversity and unique morphologies (Prangishvili et al., 2017; Munson-McGee et al., 2018; Wirth and Young, 2020). In this regard, further in vivo characterization of the putative Caldarchaeales viruses is needed to better understand their ecological roles and elucidate the co-evolution of them with their hosts in extreme habitats.

CRISPR-Cas systems in Caldarchaeales genomes have been previously reported, with MAGs harboring type-I and type-III CRISPR-Cas systems (Nunoura et al., 2011; Hua et al., 2018). Our analysis shows that type-I and type-III CRISPR-Cas systems might be also present in the MAGs reported here. We found that P. maunauluense has two cas loci classified as subtype III-D and putative subtype I-A, respectively, in two contigs (Figure 5; Supplementary Table S8). However, neither has a nearby CRISPR array (Supplementary Table S9). In the third contig of P. maunauluense, we detected two CRISPR loci, both of which lack cas genes in their vicinity. Calditenuis fumarioli and 1054_113_bin.2 have both subtype III-D and subtype I-A cas loci with nearby CRISPR arrays. 1054_108_bin.3 includes a putative subtype I-A cas locus, but this MAG lacks CRISPR arrays, possibly because it is highly incomplete and fragmented. Although we did not recover any cas locus in 1054_113_bin 10, we found that it harbors two cas3 genes in the same contig, which does not contain a CRISPR array. Furthermore, three CRISPR arrays were found in three 1054_113_bin.10 contigs.

Pelearchaeum maunauluense, C. fumarioli, and 1054_113_bin.10 encode genes for the following archaeal flagellar components: FlaB, FlaG, FlaH, FlaI, FlaJ, and FlaK (Supplementary Table S3). However, 1054_108_bin.3 only has the gene encoding FlaK. No archaeal flagellar components were identified in 1054_113_bin.2, possibly due to the low level of completeness of this genome. The presence of archaeal flagellar components implies that microbial lineages represented by P. maunauluense, C. fumarioli, and 1054_113_bin.10 may be motile under certain environmental conditions. In addition, the gene encoding for an archaeal type IV pilin was identified in P. maunauluense, C. fumarioli, 1054_108_bin.3, and 1054_113_bin.10. Archaea type IV pili are involved in biofilm formation, including surface adhesion, microcolony formation, and regulation of flagella-dependent motility (Pohlschroder and Esquivel, 2015). In this respect, biofilm formation may protect against environmental stressors such as low or high pH and toxic chemicals and promote horizontal gene transfer and syntrophy with other microorganisms (van Wolferen et al., 2018). Given the role of type IV pili in archaeal biofilm formation, organisms from which P. maunauluense, C. fumarioli, 1054_108_bin.3, and 1054_113_bin.10 were derived may be able to form biofilms. Moreover, two chemotaxis-related genes (methyl-accepting chemotaxis protein and a heme-based aerotactic transducer) were found in P. maunauluense genome. Taken together, motility genes in P. maunauluense, C. fumarioli, 1054_108_bin.3, and 1054_113_bin.10, and chemotaxis genes in P. maunauluense, provide genomic evidence that these lineages might be able to move toward favorable environments within these geothermal habitats.

The pangenomic comparison of 1054_113_bin.10 with its closest relatives from GTDB showed 1,844 genes in this pangenome, 528 of which are core, 26 which are unique, and that 19 assigned to the bin have COG categories (Figure 6A; Supplementary Table S11). Most of these gene clusters belong to “Amino acid transport and metabolism,” and “Inorganic ion transport and metabolism” COG categories that comprise ABC-type nitrate/sulfonate/bicarbonate transport system components, and ABC-type dipeptide/oligopeptide/nickel transport system components. There are also two gene clusters of the COG category “Defense mechanisms” that encode toxin-antitoxin (TA) system-associated proteins. Two reference genomes from Jinze hot spring sediment environment, China (GCA_011364605.1 and GCA_011364015.1), have 265 gene clusters that were not found in 1054_113_bin.10 or the reference genome (GCA_011373365.1) from Gongxiaoshe hot spring sediment, China. Most of the 173 gene clusters assigned to COG categories belong to “Energy production and conversion,” “Amino acid transport and metabolism,” and “Defense mechanisms.” The gene clusters of “Energy production and conversion” mainly comprise genes involved in the non-phosphorylated Entner-Doudoroff pathway, TCA cycle, and pyruvate oxidation. Gene clusters of “Defense mechanisms” mostly include genes associated with toxin-antitoxin (TA) systems, CRISPR-Cas systems, and ABC-type multidrug transport systems. On the other hand, COG category “Inorganic ion transport and metabolism” contain gene clusters encoding sulfite exporters, an ABC-type Fe³⁺-hydroxamate transport system component, ABC-type cobalamin/Fe³⁺-siderophores transport system components, a phosphate uptake regulator, and predicted Fe²⁺/Mn²⁺ transporter. There is also one gene cluster of the COG category “Mobilome: prophages, transposons” that encodes Mu-like prophage FluMu protein gp28, indicating the occurrence of past viral infections. The reference genome (GCA_011373365.1) from China’s Gongxiaoshe hot spring sediment contains 71 unique gene clusters, 34 of which are annotated; most of the annotated gene clusters belong in the COG categories “Carbohydrate transport and metabolism,” “Coenzyme transport and metabolism,” and “Post-translational modification, protein turnover, chaperones” that includes gene clusters encoding chaperonin GroEL.

The pangenomic analysis of C. fumarioli and its closest relatives from GTDB revealed 1869 gene clusters in the pangenome, with 1,253 core gene clusters, and 146 gene clusters unique to C. fumarioli, 68 of which were annotated to COG categories (Figure 6B; Supplementary Table S12). The COG categories “Carbohydrate transport and metabolism” and “Defense mechanisms” contained the largest number of unique gene clusters. The “Carbohydrate transport and metabolism” category includes ABC-type sugar transport system components and a 5-Carboxyvanillate decarboxylase (LigW) that catalyzes the conversion of 5-carboxyvanillate to vanillate in the biochemical pathway for the degradation of lignin (Vladimirova et al., 2016). The “Defense mechanisms” category mainly includes CRISPR-Cas system-associated proteins, and additionally a few toxin/antitoxin (TA) system proteins, and an ATPase component of the ABC-type multidrug transport system. Considering that two reference genomes used in the pangenome analysis of C. fumarioli were obtained from a hot spring sediment metagenome and they belong to the same species as C. fumarioli, and the proteins associated with CRISPR-Cas and TA systems are among the most common unique gene clusters of C. fumarioli, viruses infecting Archaea may be more prevalent in fumaroles than in hot springs, and the unique gene clusters might reflect adaptation strategies to the fumarole environment of organisms represented by C. fumarioli.

For example, one of the unique gene clusters in the COG category “Energy production and conversion” encodes nitrite reductase (NO-forming), part of denitrification, and an “Intracellular trafficking, secretion, and vesicular transport” category unique gene cluster that encodes VirB4 ATPase of the Type IV secretion systems, which mediate the transfer of proteins and DNA across the bacterial cell membrane and play important roles in bacterial pathogenesis and horizontal transfer of antibiotic resistance (Walldén et al., 2012). Moreover, the presence of “Inorganic ion transport and metabolism” category unique gene clusters comprising components of ABC-type Fe³⁺ transport system, ABC-type nitrate/sulfonate/bicarbonate transport system, sulfite exporter TauE/SafE/YfcA and related permeases, and enzymes associated with assimilatory sulfate and assimilatory nitrate reductions, might provide insights into the adaptations of Caldarchaeales to the physicochemical characteristics of the fumaroles. Conversely, two reference genomes from the Jinze hot spring sediment share 110 gene clusters that were also not found in C. fumarioli. Among these, 54 were annotated, and a large number belong in the COG categories “Amino acid transport and metabolism,” “Translation, ribosomal structure and biogenesis,” “Nucleotide transport and metabolism,” and “Replication, recombination and repair.” One gene cluster in the COG category, “Inorganic ion transport and metabolism,” was annotated as a predicted copper/silver-translocating P-type ATPase, which suggests potential involvement in heavy metal stress resistance. However, none of these 110 gene clusters fell into the “Defense mechanisms” category, including CRISPR-Cas and toxin/antitoxin systems.

With respect to the pangenome of P. maunauluense and its closest GTDB relatives, the total number of gene clusters was 2,209, with 953 core gene clusters, and 296 gene clusters unique to P. maunauluense is 296. Of the latter, 173 were not assigned to any COG category (Figure 6C; Supplementary Table S13). Gene clusters assigned to the COG category “Defense mechanisms” comprise the majority among 123 annotated unique gene clusters. These “Defense mechanism” gene clusters consist of proteins associated with CRISPR-Cas and toxin-antitoxin systems (TAs). TAs are small genetic modules found on bacterial mobile genetic scaffolds such as plasmids, as well as on bacterial and archaeal chromosomes (Yamaguchi et al., 2011; Yang and Walsh, 2017; Song and Wood, 2020). They typically consist of two elements: a toxin that inhibits an essential cellular process, and an antitoxin that counteracts its cognate toxin (Jurėnas et al., 2022). TAs play a critical role in the distribution and evolution of bacterial antibiotic resistance (Yang and Walsh, 2017) and promote cell survival in their native habitat (Page and Peti, 2016; Yang and Walsh, 2017). It has also been demonstrated that they are involved in multi-resistant plasmid maintenance, stress management, and biofilm formation in Bacteria (Yang and Walsh, 2017).

Song and Wood (2020) reported that a primary physiological function of chromosomally encoded TA systems in bacteria is phage inhibition, that some CRISPR-Cas system elements are derived from TA systems, and some CRISPR-Cas systems mimic TA systems by reducing host metabolism to inhibit phage propagation. Given a substantial number of annotated unique gene clusters of P. maunauluense belong in CRISPR-Cas and TA systems, viruses of Archaea might be relatively more abundant in the fumarole environment from which P. maunauluense was recovered than in the marine hydrothermal vent from which two MAGs from the same genus were obtained. Other unique gene clusters of P. maunauluense comprise COG categories “Mobilome: prophages, transposons,” which includes transposases, “Inorganic ion transport and metabolism,” which includes a high affinity permease of iron and lead ions, and superoxide dismutases implicated in oxidative stress resistance, “Nucleotide transport and metabolism,” “Carbohydrate transport and metabolism,” “Amino acid transport and metabolism,” “Coenzyme transport and metabolism,” “Lipid transport and metabolism,” “Secondary metabolites biosynthesis, transport and catabolism,” “Cell motility,” “Cell wall/membrane/envelope biogenesis,” and “Energy production and conversion” which includes tungsten-containing aldehyde:ferredoxin oxidoreductases and tetrathionate reductase subunit A, “Post-translational modification, protein turnover, chaperones” that includes predicted dithiol-disulfide isomerases which might be involved in oxidative stress tolerance (Khairnar et al., 2013), “Replication, recombination and repair” that includes DNA repair photolyase, “Signal transduction mechanisms,” “Transcription,” and “Translation, ribosomal structure and biogenesis.”

On the other hand, 113 gene clusters belong exclusively to the two MAGs from the marine hydrothermal vent and were not found in P. maunauluense. Among them, 66 were assigned to COG categories. Most of the gene clusters belong to the COG categories, “Energy production and conversion,” “Amino acid transport and metabolism,” “Lipid transport and metabolism,” and “Defense mechanisms,” with the first one containing the largest number of gene clusters. “Energy production and conversion” includes perchlorate reductase subunit alpha, which catalyzes the reduction of perchlorate to chlorite, and allows anaerobic growth on perchlorate as the electron acceptor; also included is tetrathionate reductase subunit B, a tungsten-containing aldehyde:ferredoxin oxidoreductase, and heterodisulfide reductase, subunit A (polyferredoxin)/coenzyme F420-reducing hydrogenase, delta subunit complex that may be involved in flavin-based electron bifurcation. There is also one gene cluster annotated as “Phage portal protein BeeE,” which is a member of the COG category “Mobilome: prophages, transposons,” implying past viral invasions. Interestingly, the “Defense mechanisms” category does not include CRISPR-Cas system-associated proteins, but it does include TA system proteins, an antibiotic ABC transporter ATP-binding protein, and an enamine deaminase RidA (reactive intermediate deaminase A) that protects metabolic enzymes against damage by reactive intermediates (Irons et al., 2020). Some other gene clusters that might contribute to the survival of Caldarchaeales in hydrothermal vent habitats encode manganese-containing antioxidant catalase (includes spore coat protein CotJC) in the COG category “Inorganic ion transport and metabolism,” and an activator of 2-hydroxyglutaryl-CoA dehydratase (involved in the amino acid fermentation pathway) which were assigned to “Lipid transport and metabolism,” plus a subtilisin family serine protease which could display hyperthermostable features, in the “Post-translational modification, protein turnover, chaperones” category.

We should also note that even though these pangenomic analyses could offer a perspective on organism- and niche-specific features and adaptations in Caldarchaeales, it is important to consider the incompleteness of the MAGs reported here, and of the reference genomes included in the analyses when interpreting the results. Due to the incomplete nature of many closely related MAGs in publicly available databases, we were only able to include relatively small numbers of their genomes in the pangenomic analyses. It is also important to note that the third pangenomic analysis of P. maunauluense MAG involved comparisons with genomes from at least two or more different species whereas the other two pangenomic included within-species level genomes. As such, a much larger number of core gene clusters may have been identified in the multi-species pangenomic analysis.