Genomic evolution of the class Acidithiobacillia: deep-branching Proteobacteria living in extreme acidic conditions

Abstract

Members of the genus Acidithiobacillus, now ranked within the class Acidithiobacillia, are model bacteria for the study of chemolithotrophic energy conversion under extreme conditions. Knowledge of the genomic and taxonomic diversity of Acidithiobacillia is still limited. Here, we present a systematic analysis of nearly 100 genomes from the class sampled from a wide range of habitats. Some of these genomes are new and others have been reclassified on the basis of advanced genomic analysis, thus defining 19 Acidithiobacillia lineages ranking at different taxonomic levels. This work provides the most comprehensive classification and pangenomic analysis of this deep-branching class of Proteobacteria to date. The phylogenomic framework obtained illuminates not only the evolutionary past of this lineage, but also the molecular evolution of relevant aerobic respiratory proteins, namely the cytochrome bo₃ ubiquinol oxidases.

Fig. 1. Acidithiobacillia class genus- and species-level taxons defined on the basis of whole-genome relatedness indexes based on amino acid and nucleotidic sequence identity pairwise comparisons.

A Amino acid identity values calculated with CompareM (AAI^CM; Table S3a; https://github.com/dparks1134/CompareM]) vs 16S rRNA gene sequence pairwise identities (Table S3b) showing acknowledged genus (% AAI > 70%) and species (%16 S rRNA gene > 98.7%) cutoff values [63]. B Histogram of AAI% values colored by the taxonomic rank of each pairwise comparison. C Average nucleotide identity (ANI) heatmap using ANIb and ANIm alignment algorithms threshold [47, 63]. D Digital DNA–DNA hybridization index (dDDH) box plot showing divergence between strains assigned to known and novel genomic species (Table S4) with respect to the species cutoff [48, 49], [red line: % dDDH > 70%]. The acronyms presented on the right have been subsequently used throughout the manuscript.

Fig. 2. Acidithiobacillia class consensus phylogenetic trees built using the concatenated alignment of shared genes.

A The MCC Bayesian tree was obtained using MrBayes with a concatenated alignment of 16 RPs (encompassing 2416 aa) that are present in all lineages and outgroups in the analysis. The equivalent tree obtained with all the 16 RPs of the Acidithiobacillia class is presented in Fig. S1d. Alpha- and Zetaproteobacteria form basal clades, consistent with the early branching position of these classes in the evolution of Proteobacteria [93]. B The ML tree was obtained with the PhyML package using 107 conserved single-copy proteins common to 88 genomes of the Acidithiobacillia class. Genomes of seven strains (BY-02, DMC, JYC, S10, DLC-5, GGI-221, and 21-59-9) were excluded from the phylogenomic analysis due to assembly quality issues. The alignment encompassed 24,335 aa and 10,269 parsimony informative sites. C, D Single-gene trees for RpsC (S3) and RplD (L4) constructed using NJ (10000 replicates) are shown in C and D, respectively. RPs and CPs were recovered from the genome assembly versions listed in Table S1b and detailed in Tables S5 and S6. CPs were accrued with an iterative process of filtering for genome-wise occurrence, copy number variations per genome, gene integrity, and length. Sequences aligned in this figure were reduced to simplify tree representation without altering the topology found with comprehensive alignments (Fig. S1). Branches of interest are highlighted in the tree and labeled accordingly. Posterior node support is shown as fractions. Lineages are color coded according to their taxonomic affiliation as elsewhere in the manuscript.

Fig. 3. Functional analysis of the pangenome of the Acidithiobacillia class.

A Dendrogram recreating the phylogenetic relations between acidithiobacilli. Terminal nodes, representing current species-level lineages are depicted as squares, and labeled with the acronym of the species names and ancestral nodes are represented by circles. Core genome feature characteristics are displayed to the right, and pangenome characteristics to the left. B Heatmap constructed from the genome versus function presence/absence scoring matrix. Colors reflect the dose of each gene (number of orthologs per PF identified) in each genome. Data were split between gene complement compartments (core, flexible, and exclusive) to aid visualization. Strains (per species) listed are those in Table S1.

Fig. 4. Metabolic traits of the Acidithiobacillia class species.

A Phyletic patterns of relevant energy metabolism genes (and associated gene clusters). Color coding is as labeled in the figure panel. Size of the circles represents the gene dose of representative genes in each cluster. PF variants identified are detailed in Table S7c. B Functional association networks based on gene vicinity concurrence frequencies (edges), including variants per protein family (nodes) as detailed in Table S7c. Concurrence frequency (%) of each gene pair was scored for all genomes from gene annotation tables, at selected relevant genetic contexts. Statistical concurrence information was used to color the edges as indicated in the label and were used to derive functional modules. Additional details can be found in Supplementary Information (Extended Figure Legends).

Fig. 5. Occurrence, gene dose, and genetic context of protein families encoding electron acceptor functions in Acidithiobacillia class spp.

A Phyletic pattern of terminal oxidases and relevant accessory proteins. B Gene neighborhood of terminal oxidases in representative species of the class. C Maximum likelihood phylogenetic tree of the heme A synthase CtaA genetically linked to either COX or bo₃ oxidases. D Extract of the sequence alignment of the COX2 subunit and its ortholog CyoA, showing the copper-binding motif and its variation along acidithiobacilli taxa. COX2 protein CEJ18213 is used as reference to identify key residues of the copper-binding motif.

Fig. 6. Emerging evolutionary scenario for heme copper oxidases (HCO) in the Acidithiobacillia class.

Two variant A2-type (COX) and three variant A1-type (CYO) oxidases were identified in the genomes of sequenced acidithiobacilli. Left, characteristics of extant A1-2, A1-1a/1b, and A2-type HCO CoxB and CyoA subunits with respect to the conservation of the Cu_A biding motif (number of conserved residues of the six-residues Cu_A-binding site and the number of conserved cysteines in the CXXXC moiety). The A1-1a variant gene clusters are unusual in that they group together (in the same gene cluster) the deepest branching CtaA heme A synthase-encoding genes, with genes encoding late-diverging CyoAB structural subunits. Right, evolutionary transitions inferred for the HCO terminal oxidases of the Acidithiobacillia class, based on the phyletic patterns of both COX/CYO and CtaAB-encoding genes in extant lineages of the Acidithiobacillia class. Only representative species-specific patterns are depicted for simplicity. Colored circles identify the lineages by the acronyms, and colored boxes the HCO variants present in each lineage. Fully colored boxes indicate presence and opaque boxes, absence. Boxes crossed with red lines along the dendrogram, represent probable instances of gene loss in ancestral branches of the Acidithiobacillia class tree. The dendrogram is based on CP tree depicted in Fig. 2. Conservation of the cysteins of the CXXXC Cu_A-binding motif and gene context-based association with heme A biosynthesis genes is symbolized in the right panel accordingly.