New perspectives on an old grouping: The genomic and phenotypic variability of Oxalobacter formigenes and the implications for calcium oxalate stone prevention

Abstract

Oxalobacter formigenes is a unique bacterium with the ability to metabolize oxalate as a primary carbon source. Most kidney stones in humans are composed of calcium and oxalate. Therefore, supplementation with an oxalate-degrading bacterium may reduce stone burden in patients suffering from recurrent calcium oxalate-based urolithiasis. Strains of O. formigenes are divided into two groups: group I and group II. However, the differences between strains from each group remain unclear and elucidating these distinctions will provide a better understanding of their physiology and potential clinical applications. Here, genomes from multiple O. formigenes strains underwent whole genome sequencing followed by phylogenetic and functional analyses. Genetic differences suggest that the O. formigenes taxon should be divided into an additional three species: Oxalobacter aliiformigenes sp. nov, Oxalobacter paeniformigenes sp. nov, and Oxalobacter paraformigenes sp. nov. Despite the similarities in the oxalyl-CoA gene (oxc), which is essential for oxalate degradation, these strains have multiple unique genetic features that may be potential exploited for clinical use. Further investigation into the growth of these strains in a simulated fecal environment revealed that O. aliiformigenes strains are capable of thriving within the human gut microbiota. O. aliiformigenes may be a better therapeutic candidate than current group I strains (retaining the name O. formigenes), which have been previously tested and shown to be ineffective as an oral supplement to mitigate stone disease. By performing genomic analyses and identifying these novel characteristics, Oxalobacter strains better suited to mitigation of calcium oxalate-based urolithiasis may be identified in the future.

Keywords: Oxalobacter formigenes; gut microbiome; kidney stone disease; nephrolithiasis and hyperoxaluria; oxalate degradation; phylogenomic and comparative genomic analyses; revised taxonomy.

Figure 1

Phylogenetic evaluation of the Oxalobacter formigenes taxon. (A) Maximum likelihood phylogenetic tree constructed from 691 aligned single-copy core orthologues with B. cepacia ATCC 25416 as the outgroup. Greyscale color represents original grouping and colored bars denote new species designations. Scale bar designates mean number of amino acid substitutions per site. Blue circles denote nodes with greater than 80% bootstrap support. (B) Heatmap depicting average nucleotide identity pairwise comparisons of strains used in this study. (C) Box and whisker plot displaying the whole genome percent GC content of each of the Oxalobacter species in this study. Boxes represent first and third quartile values while black line denoting the median, and whiskers encompass maximum and minimum values. (D) Coding sequence (CDS) presence/absence plot generated from pangenome. Each column represents the intersection of CDSs in that group (denoted by the number above the column). Filled circles show strains part of the same group. The bars on left show the total number of CDSs present in each genome.

Figure 2

16S rRNA gene phylogenetic analysis of the Oxalobacter genus. (A) Maximum likelihood phylogram of aligned 16S rRNA gene sequences from Oxalobacter with B. cepacia ATCC 25416 as the outgroup. Blue circles denote nodes with greater than 50% bootstrap support. Scale bar designates mean number of nucleic acid substitutions per site. (B) Alignment of consensus 16S rRNA gene sequences from each species of Oxalobacter. Red regions indicate nucleotides different than the than the consensus sequence, which excludes the V3-V4 and V4 primers. Dark gray regions show gaps in the gene. Light gray regions show similarity. Below is a nucleotide level view of the 445–472 bp region of the 16S rRNA gene.

Figure 3

Functional annotation of Oxalobacter genomes. (A) Principal coordinate analysis (PCoA) plot with Bray–Curtis dissimilarity calculated from absolute abundance of cluster of orthologous group (COG) categories. (B) Heatmap displaying that absolute abundance of each COG category. (C) Gene arrow maps of Cas clusters found in Oxalobacter species. Note cas8c gene was not identified in O. aliiformigenes OxK and the trailing cas1 gene was not identified in O. formigenes HOxHM18. (D) Circular heatmap showing the abundance of bacteriocins and antibiotic resistance genes found in the analyzed Oxalobacter genomes. (E) Heatmap showing the presence and completeness of prophages found in the analyzed Oxalobacter genomes.

Figure 4

In silico analysis of oxalyl-CoA decarboxylase. (A) Consensus alignment of the oxalyl-CoA decarboxylase amino acid sequence for each Oxalobacter species. Red regions indicate nucleotides different than the consensus sequence. Dark gray regions show gaps in the gene. Light gray regions show similarity. (B) Three dimensional in silico ribbon model reconstruction of the oxalyl-CoA decarboxylase. Also shown is the binding of the natural ligand, oxalyl-CoA, in a ribbon and space filling structure. (C) Predicted binding affinity of each model protein with the natural ligand. Data are displayed as median predicted binding affinity (kcal/mol) and analyzed by one-way ANOVA. In box plot diagrams, circles represent data point, boxes represent first and third quartile values while black lines denote medians, and whiskers encompass maximum and minimum values. ns, not significant.

Figure 5

Microbiota analysis of oxalate samples inoculated with chemostat culture. (A) Principal component analysis (PCA) plot of ASVs from sample bacterial communities. Centre log ratio-transformed Aitchison distances of sequence variants were used as input values for PCA analysis. Distance between samples on the plot represents differences in microbial community composition. Approximately 42% of the total variance is explained by the first two components shown. Strength and association for sequence variants are depicted by the length and direction of arrows shown. Points are colored by oxalate exposure and shapes denote time point (n = 3). (B) Relative abundance bar plot of longitudinal samples. Each vertical bar denotes relative sequence variant abundance (collapsed at genus-level identification). Bars are grouped by time point. (C) Analysis of differences in sequence variant abundance with oxalate supplementation controlling for the variation associated with time. Positive values indicate sequence variants that were increased with oxalate supplementation and negative values indicate sequence variants that were decreased. Taxa are defined by the most accurate level of classification. Effect size was computed with ANCOM (W ≥ 0.7) and supported with MaAsLin2 (Benjamini–Hochberg adjusted p < 0.05). (D) Longitudinal relative abundance of Oxalobacter spp. Data displayed as mean ± SD. Statistical comparisons using two-way ANOVA with Šídák’s multiple comparisons test (n = 3). (E, F) Relative abundance of predicted oxalate degrading KEGG Orthology numbers over time. Data represents mean ± SD. (E) Sum of all predicted oxalate degrading KEGG Orthology numbers (K08177, K01596 [EC:4.1.1.2], K01577 [EC:4.1.1.8], K07749 [EC:2.3.8.16], and K18702 [EC: 2.8.3.19]). (F) Abundance of K08177 (major facilitator superfamily transporter, oxalate/formate antiporter family, oxalate/formate antiporter). Statistical comparisons using two-way ANOVA with Šídák’s multiple comparisons test (n = 3). *, p < 0.05.

Figure 6

Species-specific growth of Oxalobacter in a simulated fecal environment. Oxalobacter species were individually inoculated into a simulated fecal environment and subsamples were taken every day for 5 days. (A) Relative abundance and (B) rate of change in relative abundance of each strain over time as detected by species-specific qPCR analysis. Data are displayed as mean ± SD.