To discover potentially novel Acidithiobacillia class Clade 1 B/C bacteria, we downloaded the 16S small subunit ribosomal RNA gene sequences from GenBank (as of April 2023), showing sequence identity of >97% to the reference gene sequence recovered from the ‘Igneacidithiobacillus copahuensis’ VAN18-1 strain (locus: HFQ13_RS04740). Recovered 16S rRNA gene sequences were oligotyped, according to Nuñez et al. (2017). Accession numbers and annotations for the dataset retained for further analysis are presented in Supplementary Table S1. Public genomes of Clade 1 B/C strains (VAN18-1: JAAXYO01; VAN18-2: JAAXYU01; VAN18-4: JAAXYS01; CV18-2: JAAXYP01; CV18-3: JAAXYQ01; BN09-2: JAAXYR01; YTS05: CP094359.1, CP094360.1; S30A2; JALQCS01) and reference strains used as outgroups and controls (`Fervidacidithiobacillus caldus’ ATCC 51756<sup>T</sup>: CP005986-CP005989; `Ambacidithiobacillus sulfuriphilus´ DSM 105150^T: RIZI01 and Thermithiobacillus tepidarius DSM 3134^T: AUIS010) were obtained from NCBI¹ on April 2023.

We recovered the MAGs assigned by Moya-Beltrán et al. (Moya-Beltrán et al., 2021) to ‘Ca. Igneacidithiobacillus yellowstonensis’ (Spst-908: formerly DTMS01, now SRR7540054) and ‘Ca. Igneacidithiobacillus taiwanensis’ (UBA2486: DDOU01), as well as three novel MAGs of interest identified as ‘Igneacidithiobacillus’ representatives in hot spring metagenomes of the Taupo Volcanic Zone (BioProject: PRJNA644733) with the following NCBI whole genome sequence identifiers: TVZ_G2, JAEPKW01; TVZ_G3, JAEPKX01 and TVZ_G4, JAEPKY01.

MAG sequence mCHCt20-1 (JAWNZB01), included in this study, was recovered from environmental sequencing data generated for a slurry sample collected in 2020 from a hot pool at ChanchoCó (BioBio region, Chile) designated CHCt (−37.818611 S, −71.163611 W; 1,798 m.a.s.l.; 38.9–56°C; pH 5.8–7.0). Sample manipulation for DNA sequencing and sequence data pre-processing and assembly were performed as detailed in the study by Degli Esposti et al. (2023). Contigs larger than 1,000 bp were grouped into genome bins by MaxBin2 v2.2.7 (Wu et al., 2016), Metabat2 v2 (Kang et al., 2019), and CONCOCT v1.1.0 (Alneberg et al., 2014) using the default parameters. Conserved marker genes were evaluated using CheckM v1.1.3 (Parks et al., 2015), CheckM2 v0.1.2 (Chklovski et al., 2023), and QUAST toolkit (Mikheenko et al., 2018). MAG assembly statistics were used to classify MAGs according to accepted quality standards (Bowers et al., 2018). Small-subunit rRNA sequences in contigs were identified using Barrnap (Seeman, 2017).

To evaluate the taxonomic assignment of strains and MAGs within the assembled set, genomic indexes based on amino acid and nucleotide sequences data were derived from the available genomic information. The average amino acid identity (AAI) was calculated using the CompareM program² and the aai.rb implementation from the Kostas Lab (Rodriguez-R and Konstantinidis, 2016; https://github.com/lmrodriguezr/enveomics; commit signature: fae592f) and run using the default parameters. Genus- and species-level AAI% cutoff thresholds used were < 62–70% and > 95.38, respectively (Richter and Rosselló-Móra, 2009). The average nucleotide identity based on Blast (ANIb) or Mummer (ANIm) as alignment algorithm, were calculated for all possible genome-MAGs pairs using a Python module implemented by Pritchard et al. (2016), which was available at https://github.com/widdowquinn/pyani. The in silico DNA–DNA hybridization index (dDDH) was assessed using the Genome-to-Genome Distance Calculator, and the cutoff values defined by Meier-Kolthoff et al. (2013, 2014), available at http://ggdc.dsmz.de.

Gene calling and annotation were performed using the NCBI Prokaryotic Genome Annotation Pipeline (PGAP, Tatusova et al., 2016; Haft et al., 2018). Low-quality MAGs (DDOU01 and DTMS01) were annotated through the Rapid Annotation Subsystem Technology pipeline (RAST, Aziz et al., 2008). Functional assignments were validated against the Conserved Domain Database v.3.16 (CDD, Marchler-Bauer et al., 2017) using CDsearch (Marchler-Bauer and Bryant, 2004) and hhsearch (Fidler et al., 2016) with default parameters. All predicted proteins were also analyzed against the profiles stored in the KEGG and COG databases (updated on July 2023) using the SqueezeMeta pipeline (Tamames and Puente-Sánchez, 2019). Annotated proteins were clustered by identity (threshold of 50% and default parameters) using Usearch v1.2.22 (Edgar, 2010). Protein alignments were constructed using the MAFFT v7.123 program (Katoh and Standley, 2013). Selected clusters and proteins were re-annotated and curated manually.

Gene orthology was calculated by the GET_HOMOLOGUES software package v3.3.2 (Contreras-Moreira and Vinuesa, 2013) using COGtriangles v2.1 as clustering algorithm (Kristensen et al., 2010), to recover Protein Families (PFs). BLAST pairwise alignment cutoffs were set at 75% coverage and, E-values were set at 10E-5. Gene copy numbers per PF and genome/MAG were scored by in-house Perl/Phyton scripts using E-value cutoff of 10E-5 and identity cutoff of 60%, to discriminate orthologs from paralogs. The results were curated based on genomic context analysis and manually revised alignments, when appropriate. Pan-genome metrics, including the size of core, flexible, and exclusive gene complements, were calculated by Tettelin et al. (2005). Phyletic patterns were constructed by scoring the presence/absence of representatives of each PF in each genome/MAG using parse_pangenome_matrix.pl. and other in-house scripts. Data analysis and visualization were conducted using the R package tidyverse v1.3.0.

Small subunit ribosomal RNA gene sequences and different sets of protein families of ecophysiological/evolutionary interest (Fis; conserved core proteins) were aligned using MAFFT v7.310 software with the L-INS-I method (Katoh and Standley, 2013). The 16S rRNA gene alignments were trimmed and masked (>50%) using trimAl v1.2 (Capella-Gutiérrez et al., 2009) and checked manually. Protein alignments were generated individually for each protein family with MAFFT settings maxiterate 1,000 and localpair, and the resulting alignments were trimmed using trimAl v1.2 with gap threshold of 0.5 and further refined manually. Phylogenetic trees were generated by neighbor-joining (NJ), maximum likelihood (ML), and Bayesian inference (BI) analysis to produce alternative phylogenies (Moya-Beltrán et al., 2021). The 16S rRNA gene NJ algorithm was implemented in quicktree (Howe et al., 2002), with the following settings upgma and kimura, and bootstrap values were calculated over 10,000 iterations. The ML tree was constructed using PhyML v3.0 (Guindon et al., 2010), with the following settings: General Time Reversible model (GTR, Tamura and Nei, 1993) was used as substitution model, PhyML estimated the transition/transversion ratio, the proportion of invariant nucleotides, and a discrete gamma approximation with k = 4 and 10,000 bootstrap replicates. The topology of the tree and the length of the branches were optimized by PhyML using Nearest Neighbor Interchange and Subtree Pruning and Regrafting. The Bayesian analysis tree was built with MrBayes v3.2.7 (Huelsenbeck and Ronquist, 2001) and run for 3,000,000 generations, saving trees every 100 generations. Posterior probabilities were calculated after discarding the first 30% of trees. Maximum likelihood Fis and core protein trees were reconstructed by PhyML v3.0 using Le Gascuel (LG) and Whelan and Goldman (WAG) amino acid substitution models, respectively, and tree topology optimization was performed using SMS v1.8.4 with 10,000 bootstrap replicates. Bayesian phylogenetic analysis was conducted on the core and Fis protein sets using MrBayes v3.2.7. Alignments were trimmed and manually curated before phylogenetic inference. BI trees were constructed using the Whelan and Goldman (WAG) and Jones-Taylor-Thornton (JTT) amino acid substitution models. The settings used were as follows: pRSET was used to select the amino acid substitution model, with a posteriori probability deviation equal to or near zero and lset for gamma-distributed rate variation. Each analysis was executed with MCMC to set the number of generations to 10,000, with a sample frequency of 100 generations.

Pandas (Reback et al., 2020) and Scikit-learn (Pedregosa et al., 2011) python libraries were employed for metadata-based clustering of all strains studied (see first section of Results). K-means was performed using cluster purity as the error function chosen to select the best number of clusters, altering only the number of clusters. A Random Forest algorithm was trained to predict the most common lineages. Metadata features were ranked by their ability to separate lineages. K-means was used with the top scoring features and the same purity selection protocol. In addition, the chi-square test was also tried to select top-ranked strain attributes. Both selection criteria were applied to those features available for all strains and, in a second run, to all features annotated for at least 100 strains. Visual inspection of the clustering, before and after selecting features, was performed with Seaborn clustermaps (Waskom, 2021).

Summary statistics and figures were computed using the R packages: gdata v2.18.0, dplyr v1.0.2, plotly v4.9.0, ggplot2 v3.2.1, scales v1.0, RColorBrewer v1.1.2, readr v1.2.1, and Rbase v3.6.1 implemented in Rstudio v1.2.50001 (RStudio Team, 2020). Gene cluster comparisons were obtained with clinker & clustermap (Gilchrist and Chooi, 2021). Improvement in vectorial figures was made using Inkscape v1.3.2.³ FigTree v1.4.4 was used for tree visualization and manipulation.⁴

To assess variability of the Clade 1B/C phylotypes in the sampled and sequenced space, a comprehensive phylogenetic tree was constructed with a set of 82 16S rRNA gene sequences recovered from GenBank (April, 2023; Supplementary Table S1). The genes analyzed include orthologs of the small ribosomal subunit 16S rRNA recovered from genomes and MAGs, which are tentatively assigned to ‘Igneacidithiobacillus’, when available. The neighbor-joining (NJ) phylogenetic tree built with this dataset is shown in Figure 1A. The alignment encompassed 1,272 bp of the full 16S rRNA gene sequence, with 223 variable sites and 82 parsimony informative sites (see also Supplementary Figure S1). The NJ tree obtained and rooted with the 16S rRNA gene of Thermithiobacillus tepidarius DSM 3134^T, showed a clear separation of all Clade 1 sequenced representatives from three distinctive subclades, one of which groups exclusively ‘Fervidacidithiobacillus caldus’ strains (clade 1A, represented by ATCC 51756^T). The other two sister subclades share a more recent common ancestor, each encompassing approximately half of the sequences involved in the analysis. Clade 1B includes strains originating in sulfidic caves from Mexico and sulfide ores from Colombia and China (Supplementary Table S2); in turn, Clade 1C includes strains from geothermal sites across the world (e.g., VAN18-1 from the Copahue Volcano, Argentina-Chile and V1 from Vulcano, Italy). In addition, sister subclades with diverse node depth became apparent in both branches of the tree, pertaining to ‘Igneacidithiobacillus’ (Table 1). Among these 16S rRNA sequences, 62% pertained to sequence clones from uncultured representatives of the class, 13.3% of which came from sites lacking verifiable metadata. Approximately 38% corresponded to isolates, of which only 8 strains (17.4% of cultured) had drafts or fully sequenced genomes (Supplementary Table S1). Inadvertently, some unclassified or misclassified strains of this group have been acknowledged for the biotechnological potential in applications ranging from the removal of sulfur-containing malodorous gases to chalcopyrite bioleaching (Lee et al., 2006; Kumar et al., 2008; Feng et al., 2012; Zhu et al., 2014; Li et al., 2018).

Given that certain MAGs assigned to ‘Igneacidithiobacillus’ (Moya-Beltrán et al., 2021) or inferred to pertain to the genus (TVZ_G2-G4; Sriaporn et al., 2021) lack the 16S rRNA gene marker in their current assemblies, we analyzed the tree for the fis gene product, which has recently been proven to be phylogenetically informative for the Acidithiobacillia class (Beard et al., 2023). The Fis ML tree produced with proteins recovered from the protein repository, sequenced genomes, MAGs tentatively assigned to ‘Igneacidithiobacillus’ (Supplementary Table S3), and similarly rooted (T. tepidarius DSM 3431^T), allowed us to further resolve the branching pattern (Figure 1B). In this case, Clade 1B and 1C strains appeared in sister clades encompassing an unclassified species (1B, represented by strain S30A2) and proposed ‘Igneacidithiobacillus copahuensis’ representatives (1C-I, 7 strains) and ‘Ca. Igneacidithiobacillus yellowstonensis’ single MAG (1C-II, SpSt-908). ‘Ca. Igneacidithiobacillus taiwanensis’ single MAG (UBA2468) was branched apart from the clade 1B/C together with the TVZ_G2 MAG from Tikitere in New Zealand, reported by Sriaporn et al. (2021) (novel clade 1E-I), and a MAG from ChanchoCó in Chile, mCHCt20-1 (novel clade 1E-II), contributed in this study. One additional clade, more ancestral in origin, accommodates the TVZ_G3 (novel clade 1F-I) and the TVZ_G4 (novel clade 1F-II) MAGs from Tikitere in New Zealand.

Collectively, these results indicate that additional taxa pertaining or related to ‘Igneacidithiobacillus’ exist in the currently sampled dataset, further supporting the recognition of this taxonomic unit as a genus of the Acidithiobacillia class. If the novel Fis clades (1E-I, 1E-II, 1F-I, and 1F-II) represent species of the genus ‘Igneacidithiobacillus’, or additional genus-level rank taxa, remains untested.

To evaluate if novel clades 1E-I, 1E-II, 1F-I, and 1F-II, uncovered through single marker phylogenetic analysis, represented species of the genus or different taxonomic rank units, we performed pairwise genomic comparisons between all genomes and MAGs, and derived the average amino acid identity (AAI) and the average nucleotide identity (ANIb) indexes (Figure 2; Supplementary Table S4). All of the AAI values obtained from the pairwise comparisons exceeded the recognized thresholds for genus-level differentiation (<65–72%, Konstantinidis and Tiedje, 2007), except for control taxa ‘F. caldus’ (average of all igneacidithiobacilli vs. ‘F. caldus’ = 69.2%) and T. tepidarius (average of all igneacidithiobacilli vs. T. tepidarius = 59.1%), which were lower than these threshold values (Figure 2A). These results imply that all genomes and MAGs of novel Fis-tree clades 1E-I, 1E-II, 1F-I, and 1F-II pertain to the ‘Igneacidithiobacillus’ genus along with representatives of clades 1B and 1C. In turn, the ANI data (Figure 2B) revealed seven genomic clusters below the 95.9% species threshold, all of which also emerged from the dDDH analysis using the 70% threshold (Supplementary Table S4). Therefore, genomic relatedness indexes support the existence of four novel genomic species within the ‘Igneacidithiobacillus’ genus (Supplementary Table S5), which add to the three currently acknowledged species, namely, ‘Igneacidithiobacillus copahuensis’ (typified by strain VAN18-1^T), ‘Ca. Igneacidithiobacillus yellowstonensis’ (typified by MAG SpSt-908^TS as sequence type material), and ‘Ca. Igneacidithiobacillus taiwanensis’ single MAG (typified by MAG UBA2468^TS). Proposed names for these novel species-level rank taxa are: species ‘Igneacidithiobacillus siniensis’ (typified by strain S30A2^T), ‘Ca. Igneacidithiobacillus chanchocoensis’ (typified by MAG mCHCt20-1^TS), ‘Ca. Igneacidithiobacillus taupoensis’ (typified by MAG TVZ-G3^TS), and ‘Ca. Igneacidithiobacillus waiarikiensis’ (typified by MAG TVZ-G4^TS). The multiprotein phylogenetic trees shown in Figure 2C; Supplementary Figure S2, which are constructed with a concatenate of 87 universally conserved ribosomal proteins and other phylogenetically informative housekeeping class-core marker genes, and fully supports the differentiation of type strains and type sequence material from the named sulfur-oxidizing species of the Acidithiobacillia class and acknowledged species of the genus ‘Igneacidithiobacillus’ (Supplementary Table S6).

Basic genomic features of the acknowledged species representatives are shown in Table 2. Data conveyed from strains and MAGs (average estimated completeness 87.9%) show general consistency in terms of size (average 2.4 Mb) and G + C content (average 58.1%), with most species representatives of the genus having fairly small genomes compared with the class average (3.1 Mb, Moya-Beltrán et al., 2021). I. copahuensis and I siniensis, both species with cultured representatives (average completeness 98.1%), have genome sizes reaching 2.8 Mb, which are comparable to those of its closest genus ‘Fervidacidithiobacillus’ (average size ‘F. caldus’ 2.86 Mb, Valdes et al., 2009; You et al., 2011; Moya-Beltrán et al., 2021). Larger differences and G + C content were observed for MAGs, TVZ-G3, and TVZ-G4.

To infer the ecological niches occupied by the recognized Candidate species, we evaluated the global occurrence and distribution of members of the ‘Igneacidithiobacillus’ genus. Environmental and geographical data; associated with each strain, 16S rRNA sequence clone, and/or metagenome-derived genome included in the study, was recovered from the published literature and public databases and used in statistical analysis (Supplementary Table S2). Over 95% of all igneacidithiobacilli strains and sequence clones sampled were mapped to sites along the Ring of Fire in Asiatic countries (China, Taiwan, Korea, and Japan), New Zealand, Chile, Argentina, Colombia, Mexico, and the United States (Figure 3A; Supplementary Table S2), in acknowledged tectonically or geothermally active areas (Renaut and Jones, 2011). Additional origins for ‘Igneacidithiobacillus’ spp. could be traced to volcanic sites away from the Ring of Fire such as the Vulcano island in Italy and geothermal areas in central Africa or to sulfide-rich and sulfate-rich waters and soils inland in a number of locations in China, India, and Hawaii. No obvious trend in the geographic distribution of the strains/clones assigned to the different Candidate species or phylotypes could be derived from the observed distribution map (Figure 3A) other than the more ubiquitous and cosmopolite occurrence of clade 1B and 1C-I representatives with respect to those of clades 1E and 1F (and probably attributable to sampling biases).

In the analyzed dataset, approximately half of the representative sequences originated from natural environments, such as karst caves, volcanic rivers, sinters, and geothermal pools. The other half of the ‘Igneacidithiobacillus’ spp. representatives came from anthropogenic environments (such as mines—heaps, tailings, acid mine drainages, and bioreactors—tanneries, wastewater treatment plants, and sewages) or heavily impacted environments (such as contaminated/polluted water courses or soils). The type of sample from which isolates and clones were recovered included water, sludge, sediment, soil, and snottite. All of these are recognized oligotrophic environments, which are rich in either sulfur, sulfide, or their oxidations products (Parker et al., 2008; Varekamp et al., 2009), and most are also polluted with heavy metals.

To obtain a comprehensive view of the distinctive phenotypic properties of the group, we next compared the available genomes and MAGs of ‘Igneacidithiobacillus’ spp. (n = 14, Table 2) and reconstructed relevant aspects of their metabolism. To this end, proteins were grouped by orthology into protein families, and their occurrence and conservation were evaluated against the class-wide pangenome. The igneacidithiobacilli (encompassing 14 genomes of 7 species) harbored 5,415 distinct PFs, which was significantly lower than the pangenome size of the monospecific genus ‘Fervidacidithiobacillus’ comprising 7,127 PFs across 18 genomes (Moya-Beltrán et al., 2021). While the reduced number of PFs could be attributed to the inclusion of several MAGs potentially lacking some genes, genomic completeness deduced from the occurrence of 258 universal protein markers (as implemented in CheckM, Parks et al., 2015) averaged 98.15% for genomes and varied from 66.92 to 95.37% for MAGs, with five of the seven MAGs meeting the established high-quality threshold (> 90% completeness, Bowers et al., 2018). In addition, the core/pangenome ratio of these two genera was comparable (‘Igneacidithiobacillus’: 16.5%; ‘Fervidacidithiobacillus’: 17.2%; Moya-Beltrán et al., 2021), implying that the accessory gene complement of the igneacidithiobacilli is currently under-sampled.

A total of 896 PFs were shared (core) among representatives of the seven ‘Igneacidithiobacillus’ lineages (Figure 4A). This core represents 31.6% of the predicted proteome of `I. copahuensis´ VAN18-1 type strain (2,838 CDSs; 2,749 PFs), being slightly higher than the percentage of core proteins carried by other Acidithiobacillia (e.g., 21.4% in Acidithiobacillus spp.; Moya-Beltrán et al., 2021). The number of shared genes rose in correlation with increasing genetic relatedness among the species, amounting to 1,043 for clade 1B-1E, 1,245 for clade 1B-1C members, and 1,615 for `I. copahuensis´ strains (clade 1C-I), the species of the genus with the largest number of cultured and sequenced representatives (Figure 4A). These numbers reflect relevant gene gains and losses that occurred during the differentiation of these lineages. Core genes partitioned into 21 COG categories and 25 KEGG paths and included most housekeeping genes required for basal metabolism of Acidithiobacillia class members (Moya-Beltrán et al., 2021), encompassing essential information processing genes, chemolithoautotrophic energy metabolism, and cell envelope biogenesis and maintenance genes (Figure 4B; Supplementary Table S7).

Sulfur metabolism enzyme-encoding genes required for the utilization of hydrogen sulfide (sqr), elemental sulfur (sdo, sor), and thiosulfate (sox) as electron donors, as well as those encoding key respiratory chain elements (cydA-CIO, cyoB, nuo-complex, and atp-complex), were present and conserved with respect to known acidithiobacilli (Figure 4C), strongly suggesting that the igneacidithiobacilli are also obligate chemolithotrophs. Some of these genes or complexes varied in copy number, with enrichment of functions in particular lineages (e.g., 4–5 copies of sdo in ‘I. copahuensis’) and depletion of functions in others (e.g., 1 copy of the nuo complex in ‘Ca. I. taupoensis’ and ‘Ca. I. waiarikiensis’).

As other Acidithiobacillia class members, ‘Igneacidithiobacillus’ seem to be incapable of dissimilatory sulfate reduction, lacking key functional markers for this process (dsrABC, Thauer et al., 2007). They also lack the sreABCD genes that enable anaerobic growth on sulfur coupled to ferric iron reduction in iron/sulfur-oxidizing Acidithiobacillus spp. (e.g., Osorio et al., 2013). Yet the MAGs of lineages ‘Ca. I. chanchocoensis’ and ‘Ca. I. waiarikiensis’, encode a thiosulfate/polysulfide reductase (phs/psr) known to catalyze the stepwise reduction of elemental sulfur and the zero-valent sulfur of thiosulfate or polysulfides to sulfide in the absence of oxygen (Heinzinger et al., 1995; Jormakka et al., 2008).

Certain igneacidithiobacilli (‘Ca. I. yellowstonensis’, ‘I. siniensis’, ‘Ca. I. taupoensis’, and ‘Ca. I. waiarikiensis’) seem also unable to derive energy from tetrathionate, lacking the genes encoding not only the tetrathionate hydrolase (tetH) but also the membrane-bound tetrathionate forming thiosulfate:quinone oxidoreductase (doxDA). Lineages missing these PFs are the ones originating from the highest temperature habitats, where tetrathionate would be particularly unstable in the presence of strong reductants such as sulfide (Xu et al., 1998), and which is highly abundant in most of these habitats. Genes for the assimilation of sulfur into the cysteine biosynthetic pathway identified in the ‘Igneacidithiobacilllus’ spp. included those involved in the acquisition of sulfur from sulfate (sat, cysJ, cysE, and cysM) as other Acidithiobacillia class species (Valdés et al., 2003), and from sulfonates (ssuD); so far, the capacity has been overlooked in this group of acidophiles (Supplementary Table S7). Lack of the APS kinase (cysC) and the sulfite reductase (cysJI) suggests that thiosulfate—instead of sulfate—may be used as sulfide donor in the synthesis of cysteine by CysM (Zhao et al., 2006).

Several core PFs (n = 29) occurred in gene dosages are higher than 2, including relevant functions for carbohydrate transport and metabolism (n = 8) involved in glycolysis (pdhB, pdhC, pgi, and gpmI), pyruvate metabolism (ackA), pentose phosphate metabolism (zwf), carbon fixation through the Calvin-Benson-Bassham reductive pentose phosphate cycle (rbcL), glycogen utilization (glgP, cga), and energy metabolism (cydAB), among others (Supplementary Table S7, core PFs). Several of these genes provide the igneacidithiobacilli with the capacity to fix carbon dioxide autotrophically.

Exclusive gene complements (present in all genomes/MAGs) of each lineage were predominantly hypothetical functions and variants of PFs, existing in other lineages, providing little insight into the adaptive features, and differentiating these species. Lineage-specific PFs with predicted functional assignments included genes involved in: (a) defense against foreign DNA which were particularly overrepresented in ‘Ca. I. yellowstonensis’, ‘I. siniensis’, and ‘Ca. I. chanchocoensis’, (b) osmotolerance preservation which potentiated different physiological strategies in each species, (c) cell envelope integrity and modification which entailed a wide diversity of species-specific glycosyltransferases, and (d) both uptake and efflux transporters of different sorts in each lineage.

To ascertain the physiological characteristics of the novel ‘Igneacidithiobacillus’ spp., we analyzed common and differential genome-derived traits in the light of emerging trends in the published metadata for the reassigned strains and clones (Supplementary Figure S3). Preferential temperature, pH, and water chemistry of the selected representatives of the igneacidithiobacilli lineages derived from available reports in the literature (and/or generated herein) on the source habitat, isolation conditions, or growth experiments are shown in Table 3. Clade 1C representatives (‘I. copahuensis’, ‘Ca. I. yellowstonensis’) have been categorized as thermotolerant (Satoh et al., 2006; Hassan et al., 2010; Li et al., 2018; Moya-Beltrán et al., 2021), while clade 1B representatives (‘I. siniensis’) have been reported to endure extremely low pHs (Ni et al., 2008; Feng et al., 2012; Jones et al., 2016; Liang et al., 2019). In turn, clade 1E lineages originate from mildly acidic habitats with higher temperature maxima (Lin et al., 2015), and clade 1F representatives originate from mildly thermal environments of higher pH (Sriaporn et al., 2020, 2021). The habitats of origin differ also in their reported water chemistry, having relevant differences in the contents of chloride salts, carbonates, and metals. We thus assessed the distribution of genes encoding products that could confer the igneacidithiobacilli lineages with differential growth capacities based on these habitat and growth characteristics.