Genome assemblies for currently classified Thermoanaerobacterales (based on NCBI taxonomy, March 2023) were obtained from the National Center for Biotechnology Information (NCBI). The list of species strains and accession numbers are provided in Supplementary Table S1. 16S ribosomal RNA gene sequences were retrieved based on NCBI annotations for the most complete 16S rRNA sequence found in these assemblies.

Isolation location coordinates for the Caldicellulosiruptor were obtained primarily through isolation literature and NCBI biosample information. When needed, locations were estimated as ‘best possible’ from information contained in the literature (e.g., town names or geographical features). Supplementary Tables S2, S3 contain detailed information on the locations and their accuracy. Geographical distances between locations were calculated via a webserver employing the formula developed by Vincenty (1975).¹

16S rRNA gene sequences were aligned using Clustal Omega (v. 1.2.4; Sievers et al., 2011) using flags “--outfmt fa --distmat-out --full --full-iter --percent-id --guidetree-out.” The percent identities output was used for 16S rRNA identity matrices (Supplementary Table S8), which were color coded based on proposed 16S rRNA taxonomic rank “Minimum Sequence ID” values (Yarza et al., 2014). The 16S rRNA alignment was inputted into FastTree (v. 2.1.11; Price et al., 2010) to generate a distance tree using flags “-nt -gtr -gamma.” The 16S rRNA tree was mid-point rooted with Dendroscope 3 (Huson and Scornavacca, 2012).

The bac120 gene marker set from the Genome Taxonomy Database (GTDB) was also used to infer a phylogenetic tree from the same set of Thermoanaerobacterales (based on NCBI taxonomy, March 2023) using GTDB-Tk v2.1.0 (Parks et al., 2018, 2020, 2021; Chaumeil et al., 2019, 2022). The phylogenetic tree was inferred using the “de_novo_wf” workflow with 200 aa per gene marker and using the phylum Thermodesulfobiota as the out-group for tree rooting. The GTDB reference database was excluded to allow for a custom taxonomic classification. Both the 16S rRNA and bac120 trees were visualized and formatted using the Interactive Tree of Life (iTOL; Letunic and Bork, 2006). A color-blind friendly divergent color palette was applied to the different genera (determined by Chroma.js Color Palette Helper)² (colors: #1d1b73, #243694, #2b53af, #2d74be, #009ba7, #b1ab14, #ba971c, #c5801d, #d36219, #f00000). The same genera colors were applied to the 16S rRNA tree.

Average nucleotide identities (ANI) were generated with pyani (v. 0.2.12; Pritchard et al., 2016) average_nucleotide_identity.py BLAST+ method (ANIb; flag “-m ANIb”). Final bidirectional ANIb values were generated by averaging the ANIb values representing the two orientations of the same strain pairs (Supplementary Table S6). Average amino acid identities (AAI) were calculated with GET_HOMOLOGUES (v22082022; Contreras-Moreira and Vinuesa, 2013; Vinuesa and Contreras-Moreira, 2015) with flags “-A –t 0 –M (Supplementary Table S7).” A composite ANIb/AAI matrix of the 91 genomes was calculated and heat-mapped based on the min/max values (40.4–100%) using the Inferno colormap. This same colormap scale was maintained for all ANIb/AAI figures.

Core and pangenome analysis were also completed using GET_HOMOLOGUES. For this, the flags “-A –t 0 –c –z –P –M” were used. Then compare_clusters.pl. was used to generate the core, soft core, and pangenomes. Finally, parse_pangenome_matrix.pl. with “-f core_Tettelin” was used to generate pan/core genome graphs (Supplementary Data S2).

For annotation of carbohydrate active enzymes, dbCAN (v 3.0.7; Zhang et al., 2018; Zheng et al., 2023) was locally run with the run_dbcan command using protein.faa files of translated nucleotide coding sequences as input and flags “--signalP = true --gram P” to use SignalP v 4.1 (Nielsen, 2017; Supplementary Data S2). The dbCAN output was then curated where all glycosyl transferases (GTs) were removed, and a consensus was generated between the three prediction algorithms (HMMER, eCAMI, and DIAMOND), where annotations that had agreement between at least 2 prediction algorithms was used as the consensus. Any proteins that had positive hits on only a single algorithm were subjected to a BLAST search on the NCBI webserver against the NCBI non-redundant protein sequences database to obtain an annotation. Any coding sequence deemed to not be involved in carbohydrate catabolism, based on the BLAST search annotation, were removed from the analysis, along with any coding sequences corresponding to proteins under 100 amino acids (Supplementary Table S9). For the Caldicellulosiruptorales, remaining proteins representing the “CAZysome” of each strain were input into GET_HOMOLOGUES for Pan/Core/Soft-core analysis, as described above for whole genomes. The resulting Core/Soft-cores of the Caldicellulosiruptorales were then analyzed for CAZyme annotation and predicted function (Supplementary Table S4).

Three strains of Caldicellulosiruptorales were obtained from DSMZ-German Collection of Microorganisms and Cell Cultures GmbH: DSM 6725 (C. bescii), DSM 18901 (C. hydrothermalis), and DSM 112098 (C. diazotrophicus). Strains were adapted to substrates and grown in 125 mL serum bottles, as described previously (Bing et al., 2023a) where modified D671 media with cellobiose (Biosynth-Carbosynth, OC04040) was used for routine growth. All cultures were grown at 75°C with 150 rpm shaking in a New Brunswick Innova 42 incubator shaker. Adaption to beechwood xylan (Biosynth-Carbosynth, YX45751), wet-milled corn fiber (WMCF, provided by Novozymes A/S), or Avicel PH-101 (Millipore-Sigma) was done by passages of 5 × 10⁸ cells (final 1 × 10⁷ cells/ml starting cell density) from 5 g/L cellobiose cultures to 0.5 g/L cellobiose and 4.5 g/L carbohydrate equivalent of substrate (xylan, WMCF, or Avicel). Cell growth was monitored by epifluorescence microscopy, as previously described (Bing et al., 2023a). Once cells reached 5 × 10⁸ cells/ml or after 3 days, 5 × 10⁸ cells were passaged to 5 g/L carbohydrate equivalent of substrate alone. Culture growth was monitored for 7 days. C. bescii (DSM 6725) and C. hydrothermalis (DSM 18901) were included as controls for expected growth phenotypes on the substrates. All strains were expected to grow on cellobiose, xylan, and WMCF, but only C. bescii grows on microcrystalline cellulose (Avicel). The phenotype observed is robust growth of C. bescii on all substrate adaptions and final passages; C. hydrothermalis had robust growth on cultures containing cellobiose, xylan, and WMCF, weak growth on the adaption passage to Avicel (containing 0.5 g/L cellobiose) and no growth on Avicel alone.

In order to assess the existing taxonomic classification for the genus Caldicellulosiruptor, phylogenetic trees were inferred for a set of selected genomes within the order Thermoanaerobacterales (according to NCBI Taxonomy as of March 2023) using both 16S ribosomal RNA gene sequences to represent the existing taxonomic structure and the bac120 gene markers from the Genome Taxonomy Database (GTDB) to assess reclassification (Figure 1). The larger bac120 gene marker set allowed a multiple sequence alignment of ~24,000 amino acids, enabling more robust placement of the selected genomes within the tree and highlights the several areas needed for reclassification. Note that recent publications have already moved Thermodesulfobium to its own phylum Thermodesulfobiota (Frolov et al., 2023) and the genera Calderihabitans, Desulfitibacter, Moorella, and Zhaonella to the novel order Moorellales (Lv et al., 2020); however, these changes were not yet reflected in the NCBI taxonomy at the time of this analysis.

The 16S ribosomal RNA gene sequence tree includes multiple nodes at the genus level and higher with bootstrap values ≤0.50, indicating low confidence some of the taxonomic arrangements in this tree. Meanwhile, the bac120 gene marker tree has ≥0.50 bootstrap values (and ≥ 0.98 with the exception of the Calorimonas adulescens node) for all nodes above the species level, implying strong phylogenetic relationships above species-level classification. Based on the bac120 gene marker phylogeny, the Order Thermoanaerobacterales should consist of the Family Thermoanaaerobacteraceae containing the genera Thermoanaerobacter, Thermoanaerobacterium, and Caldanaerobacter, the Family Caldanaerobieaceae containing the genus Caldanaerobius, and the Family Calorimonaceae containing the genus Calorimonas. The current genus Caldicellulosiruptor is sufficiently divergent from the order Thermoanaerobacterales that it should be long to a separate Order, Caldicellulosiruptorales, and Family, Caldicelluosiruptoraceae. Within the Family, the species C. danielli, C. hydrothermalis, C. diazotrophicus, C. owensensis, C. obsidiansis, C. bescii, C. kronotskyensis, C. acetigenus, C. lactoaceticus, and C. kristjanssonii should be moved to a new genus; we propose this genus be named “Anaerocellum” in view of the fact that the type strain C. bescii was originally named Anaerocellum thermophilum (Svetlichnyi et al., 1990; Kataeva et al., 2009; Yang et al., 2010). The species C. changbaiensis, C. sp. F32, C. saccharoloyticus, C. naganoensis, and C. morganii should remain in the genus Caldicellulosiruptor, recognizing that the type strain from this group was C. saccharolyticus (Rainey et al., 1994).

Together, the orders Thermoanaerobacterales and Caldicellulosiruptorales likely comprise their own Class (Thermoanaerobacteria), but without comprehensive evaluation of the Class Clostridia, this could not be fully established here (Supplementary Table S5). The placement of Mahella australiensis in the bac120 gene marker tree further supports the need for a broader evaluation of the Class Clostridia. The GTDB taxonomy (~5,000 amino acid alignment) keeps the genus Mahella within Clostridia, but our ~24,000 amino acid alignment suggests that Mahella is more closely related to the Caldicellulosiruporales, meriting a Family-level classification within the Order Caldicellulosiruptorales.

Reclassification at the species level is further supported by Average Nucleotide Identity (ANI) and Average Amino Acid Identity (AAI) comparisons of genome-sequenced strains currently assigned to the Order Thermoanaerobacterales (Figure 2), including 15 strains within the genus Caldicellulosiruptor (Table 1). Color coding in Figure 2B reflects previously proposed taxonomic grouping thresholds for 16S rRNA gene sequence identities (Yarza et al., 2014). Although the 16S rRNA gene sequence analysis is less clear than ANI/AAI, as expected, only a few current Thermoanaerobacterales genera are indicated as having Order level or lower relatedness by 16S rRNA.

To further investigate the proposed re-classification of the Thermoanaerobacterales, a genome-wide assessment was done (Figure 3). To investigate species-level classification, we assume that ANI values ≥94–96% are indicative of strains of the same species (Richter and Rossello-Mora, 2009; Kim et al., 2014). ANI/AAI results indicate several species from the genera Thermoanaerobacter, Moorella, Carboxythermus, Thermoanaerobacterium, Caldicellulosiruptor, and the new genus Anaerocellum need to be consolidated as strains of previously described species (Table 2). For example, Thermoanaerobacter ethanolicus (type strain JW200, the first isolate from this group) should now include an additional 6 strains currently assigned as T. thermohydrosulfuricus WC1, T. thermohydrosulfuricus DSM569, T. sp. RKWS2, T. siderophilus SR4, T. weigelii RT8.B1, and T. indiensis BSB-33, all of which share ANIb ≥97%.

Within the Caldicellulosiruptorales, several species level reclassifications are appropriate. The recently proposed reclassifying of A. lactoaceticus and A. kristjanssonii as subspecies of A. acetigenus (Habib et al., 2021) is supported by ANI values all ≥97% (Figure 3). Note that A. acetigenus was previously reclassified from Thermoanaerobium acetigenum based on 16S rRNA gene sequence and physiological properties (Onyenwoke et al., 2006). What we now term A. bescii and A. kronotskyensis should be considered the same species, given that their ANI is 96% (Table 2; Figure 3); as such, A. kronotskyensis 2002 is designated as A. bescii strain 2002. Two other kronotskyensis strains (2006, 2,902) were classified as this specices based on 16S rRNA (Miroshnichenko et al., 2008) but have no available genome sequences. These may also be strains of A. bescii, but without genome sequences this cannot be verified.

Note that the inter-genus ANI of the Caldicellulosiruptorales is ≤82%, indicative of the divergence between the genera Anaerocellum and Caldicellulosiruptor, established by the GTDB-tk analysis. The Caldicellulosiruptor genus share an ANIb >83% (89% if C. morganii is excluded; Figure 3). Note that C. morganii is the most divergent of the Caldicellulosiruptorales by ANI/AAI, although the bac120 (GTDB-tk) phylogenetic tree places it squarely within the Caldicellulosiruptor genus (Figure 1). A case can be made that C. saccharoloyticus, C. sp. F32, and C. changbaiensis CBZ are the same species (ANI’s ≥ 95% and agreement with bac120 phylogenetic inference; Table 2); in fact, C. saccharolyticus and C. sp. F32 have been designated as such by the China General Microbiological Culture Collection (CGMCC). However, the C. sp. F32 genome is currently in 127 contigs, thus lowering the confidence in its taxonomic placement compared to other Caldicellulosiruptor strains.

The core genomes of various taxonomic levels were evaluated in the process of whole-genome analysis. The core and pan genomes of the Thermoanaerobacteraceae are 306 / 6,915 for the family, 1,308 / 3,510 for Caldanaerobacter, 815 / 4,087 for Thermoanaerobacter, and 1,374 / 4,419 for Thermoanaerobacterium. Within the Caldicellulosiruptoraceae, the core and pan genomes are: 1,248 / 3,833 for the family, 1,496 / 3,527 for Anaerocellum, and 1,367 / 3,027 for Caldicellulosiruptor (Supplementary Figure S1). The pangenome of the Caldicellulosiruptoraceae remains open, implying more genetic diversity within the Family remains to be discovered.

The core and pan genomes of the Thermoanaerobacteraceae, respectively are: 306/6,915 for the family, 1,308/3,510 for Caldanaerobacter, 815/4,087 for Thermoanaerobacter, and 1,374/4,419 for Thermoanaerobacterium (Supplementary Figure S1).

The wide global distribution of the Caldicellulosiruptorales is evident from the isolation sites of currently named species as well as from signatures detected in community analyses (Figure 4; Blumer-Schuette, 2020). Given the closely related microbiological features of members of the Caldicellulosiruptorales, it is interesting to consider how these species became globally distributed and the relationships between geography, physiochemical features of isolation sites, and strain relatedness. The fact that all known members of the Caldicellulosiruptorales grow best at temperatures above 70°C differentiates them from almost all other characterized bacteria, which are mostly mesophilic or moderately thermophilic (T_opt ≤ 65°C). Community analyses of terrestrial hot springs indicate that, above 65°C, Caldicellulosiruptorales dominate (Vishnivetskaya et al., 2015; Lee et al., 2018). Based on the global presence of Caldicellulosiruptorales in low-salinity terrestrial thermal sites (>65°C) and the relatedness of strains isolated in various regions, additional areas likely to harbor novel strains of Caldicellulosiruptorales can be identified. These areas are detailed in Figure 4, both circled in the map (Figure 4A) with details provided below the map. There is a notable absence of isolated Caldicellulosiruptorales from Africa and South America, although there are thermal features on these continents hospitable to Caldicellulosiruptorales. Presence of widespread thermal features and high similarity of strains currently isolated from North America and Iceland may indicate that additional more divergent Caldicellulosiruptorales exist in these areas, which would reflect the diversity seen in China, Kamchatka, Japan, and New Zealand.

It is also interesting to consider how the globally distributed Caldicellulosiruptorales compare with respect to geographic separation of isolation sites and overall genetic relatedness. Figure 4B shows how geographical distance between isolation sites for Caldicellulosiruptorales (Supplementary Table S2) relates to AAI and ANI. While it is not surprising that species and strains isolated from immediate proximity to each other can have high AAI/ANI values, it is interesting that some species from isolation sites separated by thousands of kilometers are also closely related. Further, several strains isolated from relatively close sites have lower AAI/ANI values, indicating close proximity of isolation does not always infer high strain similarity even within the same genus. Thus, from available genomic data, there is no correlation between geographical spacing of isolation sites and ANI/AAI for the Caldicellulosiruptorales, at least at for distance >50 km. Presumably, at some smaller distance (i.e., <50 km), there could be a correlation, where multiple strains from a single geothermal feature, or very near-by thermal features, could have high relatedness. The Icelandic strains (A. acetigenus) may hint at this possibility due to their high relatedness and close proximity of isolation. However, at the same time, species isolated <100 km from each other have much lower relatedness, like those from New Zealand (C. morganii, C. saccharolyticus, A. danielii) or Japan (C. naganoensis, A. diazotrophicus).

Comparing and contrasting the underlying metabolism, physiology, and ecology of microorganisms goes beyond 16S rRNA, GTDB-tk, and ANI/AAI analyses. Interest in the Caldicellulosiruptorales was initially driven by their ability to degrade lignocellulosic biomass and the inventory of carbohydrate active enzymes (CAZymes) supporting this characteristic (Blumer-Schuette et al., 2014). Specifically, within the group of CAZymes are intracellular, extracellular, and surface-layer Glycoside Hydrolases (GHs) that process carbohydrates; many of these enzymes have multiple domains and associated carbohydrate binding modules (CBMs; Blumer-Schuette et al., 2010; Conway et al., 2016, 2017, 2018; Crosby et al., 2022; Laemthong et al., 2022). Figure 5 summarizes GH inventories separated by GH family (Drula et al., 2022), for the 15 genome-sequenced Caldicellulosiruptorales and 37 sequenced strains within the revised Thermoanaerobacterales. It is evident that the GH inventory varies widely across genera and species, and even across strains within a species. To this point, A. acetigenus str. Acetigenus has 15 or more total CAZymes, and 20 or more GH containing CAZymes, than either of the other two strains of this species.

A wide variety of GH family domains are found in all Caldicellulosiruptorales genomes; collectively these enable degradation of carbohydrates found in lignocellulosic biomass to recruit them as carbon and energy sources. Enzyme functions likely associated with these domains include xylan hydrolysis, exo- and endo-acting cleavage of β-glycans, hydrolysis of α-glucans (such as starch, pullulan), pectin degradation, hydrolysis of galactomannans, mono-, di-, and oligosaccharide phosphorylation, and hydrolysis of peptidoglycan or chitin. Domains related to cleavage of terminal glucose residues from β-glycosides (GH1) and α-glucans (GH15), as well as GH51 (likely endo acting β-glucanase or α-L-arabinofuranosidase), are found in all but one species. Domains related to phosphorylases (GH65), uncapping of glucuronic acid decorated xylooligosaccharides (GH67), and hydrolysis of unsaturated glucuronyl/galacturonyl linkages (GH105), and GH4 (possible α-glucosidase, galactosidase, glucuronidase, or galacturonase) are present in all but two species. GH9 and GH48 domains (and CBM3) are central to enzymes that catalyze microcrystalline cellulose hydrolysis, notably present in the multi-domain cellulase, CelA (Conway et al., 2017) one or both of these GH domains are present in all five Caldicellulosiruptor species and six of ten Anaerocellum species. Strains lacking these cellulases are unable to grow on microcrystalline cellulose (Avicel). Of note here, A. diazotrophicus (which lacks these enzymes) was reported in its isolation paper to grow on Avicel (Chen et al., 2021b). However, as part of this work, A. diazotrophicus was readily cultured on cellobiose, beechwood xylan, and corn fiber, but no growth was observed on Avicel at 75°C. This result reinforces that Caldicellulosiruptorales lacking GH9 and GH48 domains cannot grow on microcrystalline cellulose. Both strongly cellulolytic and non-cellulolytic (or weakly cellulolytic) species are found in most Caldicellulosiruptorales isolation biotopes (Vishnivetskaya et al., 2015; Lee et al., 2018). This reflects the cooperative nature of microbial communities. For example, not all species in a biotope need to be cellulolytic to utilize lignocellulose, but instead some species can degrade cellulose through extracellular GHs thereby enabling scavenging of the resulting oligosaccharides by non-cellulolytic species. Interestingly, common to all members of the Caldicellulosiruptorales are peptidoglycan lyases / chitinase (GH23) and peptidoglycan/ chitin binding domains (CBM50). This implicates degradation of bacterial or fungal cell wall remnants to support acquisition of substrates in otherwise nutritionally sparse biotopes.

In the Thermoanaerobacterales genomes, highly represented GH families (GH1-5, GH13, GH15, GH23, GH31, GH36, GH65, GH130) mirror what is seen in the Caldicellulosiruptorales, with two exceptions. GH18 domains (chitinases, lysozymes) are common to the Thermoanaerobacterales, but absent in the Caldicellulosiruptorales. None of the genome sequenced Thermoanaerobacterales are known to degrade crystalline cellulose, and accordingly, none contain genes encoding GH9, GH48, and CBM3 domains, which are important for microcrystalline cellulose degradation. It is also interesting that many Thermoanaerobacterales lack GH10/11 domains that relate to xylan hydrolysis (Thermoanaerobacterium, Caldanaerobius polysaccharolyticus, and some Thermoanaerobacter are exceptions). Among the Thermoanaerobacterales, Thermoanaerobacterium have the most diversity in CAZymes, reflecting their hemicellulolytic activity. Similar trends are seen with other CAZyme domains (carbohydrate binding modules, polysaccharide lysases, and carbohydrate esterases; Supplementary Figure S2). Based on CAZyme inventory, Caldicellulosiruptorales have a broader carbohydrate appetite than the Thermoanaerobacterales, especially with respect to crystalline cellulose.