Proposal of a Taxonomic Nomenclature for the Bacillus cereus Group Which Reconciles Genomic Definitions of Bacterial Species with Clinical and Industrial Phenotypes

Abstract

The Bacillus cereus group comprises numerous closely related species, including bioterrorism agent B. anthracis, foodborne pathogen B. cereus, and biopesticide B. thuringiensis Differentiating organisms capable of causing illness or death from those used in industry is essential for risk assessment and outbreak preparedness. However, current species definitions facilitate species-phenotype incongruences, particularly when horizontally acquired genes are responsible for a phenotype. Using all publicly available B. cereus group genomes (n = 2,231), we show that current species definitions lead to overlapping genomospecies clusters, in which 66.2% of genomes belong to multiple genomospecies at a conventional 95 average nucleotide identity (ANI) genomospecies threshold. A genomospecies threshold of ≈92.5 ANI is shown to reflect a natural gap in genome similarity for the B. cereus group, and medoid genomes identified at this threshold are shown to yield resolvable genomospecies clusters with minimal overlap (six of 2,231 genomes assigned to multiple genomospecies; 0.269%). We thus propose a nomenclatural framework for the B. cereus group which accounts for (i) genomospecies using resolvable genomospecies clusters obtained at ≈92.5 ANI, (ii) established lineages of medical importance using a formal collection of subspecies names, and (iii) heterogeneity of clinically and industrially important phenotypes using a formalized and extended collection of biovar terms. We anticipate that the proposed nomenclature will remain interpretable to clinicians, without sacrificing genomic species definitions, which can in turn aid in pathogen surveillance; early detection of emerging, high-risk genotypes; and outbreak preparedness.IMPORTANCE Historical species definitions for many prokaryotes, including pathogens, have relied on phenotypic characteristics that are inconsistent with genome evolution. This scenario forces microbiologists and clinicians to face a tradeoff between taxonomic rigor and clinical interpretability. Using the Bacillus cereus group as a model, a conceptual framework for the taxonomic delineation of prokaryotes which reconciles genomic definitions of species with clinically and industrially relevant phenotypes is presented. The nomenclatural framework outlined here serves as a model for genomics-based bacterial taxonomy that moves beyond arbitrarily set genomospecies thresholds while maintaining congruence with phenotypes and historically important species names.

Keywords: Bacillus anthracis; Bacillus cereus; Bacillus cereus group; Bacillus thuringiensis; bioterrorism; foodborne pathogens; phylogenetic analysis; taxonomy.

FIG 1

Dendrogram constructed using symmetric pairwise average nucleotide identity (ANI) dissimilarities calculated between 2,218 B. cereus group genomes from NCBI’s RefSeq database with N₅₀ of >20 kbp (i.e., $D_{\text{ANI}}^{\text{sym}}$ in Materials and Methods) and the average linkage hierarchical clustering method implemented in the hclust function in R. Blue tip labels denote the location of species type strain/reference genomes in the dendrogram, while tree height corresponds to ANI dissimilarity. Branch colors correspond to branch height within the tree. Dashed vertical lines appear at dissimilarities of 7.5, 6, 5, and 4, which correspond to ANI thresholds of 92.5, 94, 95, and 96, respectively (from left to right in order of appearance along the x axis).

FIG 2

Weighted undirected graphs constructed using symmetric pairwise average nucleotide identity (ANI) values calculated between 2,218 B. cereus group genomes from NCBI’s RefSeq database with N₅₀ of >20 kbp (i.e., $S_{\text{ANI}}^{\text{sym}}$ in Materials and Methods). Nodes represent individual genomes, while weighted edges connect each pair of genomes with a mean ANI value of ≥95 (A) and ≥92.5 (B), where edge weight corresponds to the mean ANI value of the pair. Nodes (i.e., genomes) are colored by (i) closest matching type strain genome or (ii) closest matching medoid genome of clusters formed at the respective ANI value. Graphs were constructed using the graphout layout algorithm implemented in R’s igraph package, using 1 million iterations and a charge of 0.02.

FIG 3

Maximum likelihood phylogenies of 2,218 B. cereus group genomes with N₅₀ of >20 kbp. Tip and branch labels are colored by genomospecies assignment using medoid genomes of genomospecies clusters formed at the widely used genomospecies threshold of 95 ANI (clusters are arbitrarily numbered) (A) and presence (colored) and absence (gray) of anthrax toxin genes cya, lef, and pagA (B); cereulide synthetase-encoding cesABCD (C); and one or more previously described Cry or Cyt insecticidal toxin-encoding genes (D). Phylogenies were constructed using core SNPs identified in 79 single-copy orthologous gene clusters present in 2,231 B. cereus group genomes. The type strain of “B. manliponensis” (i.e., the most distantly related member of the group) was treated as an outgroup on which each phylogeny was rooted. Virulence genes (cya, lef, and pagA and cesABCD) were detected using BTyper version 2.3.2 (default thresholds), while insecticidal toxin-encoding genes were detected using BtToxin_scanner version 1.0 (default settings; presence and absence of high-confidence, previously known Cry- and Cyt-encoding genes are shown, with predicted putative novel insecticidal toxin-encoding genes excluded).

FIG 4

Histogram of pairwise average nucleotide identity (ANI) values calculated between 2,231 B. cereus group genomes downloaded from NCBI’s RefSeq database. FastANI version 1.0 was used to calculate all pairwise ANI values. For histograms colored according to closest species type strain/reference genome at a conventional ≥95 ANI threshold, or histograms showing pairwise ANI values calculated between genomes meeting additional quality thresholds, see Fig. S1 and S2, respectively.

FIG 5

Maximum likelihood phylogeny of 2,218 B. cereus group genomes with N₅₀ of >20 kb. Tip and branch labels are colored by genomospecies assignment using medoid genomes of genomospecies clusters formed at the proposed genomospecies threshold of 92.5 ANI. Phylogeny was constructed using core SNPs identified in 79 single-copy orthologous gene clusters present in 2,231 B. cereus group genomes. The type strain of “B. manliponensis” (i.e., the most distantly related member of the group) was treated as an outgroup on which the phylogeny was rooted.

FIG 6

Taxonomic hierarchy for the proposed B. cereus group nomenclature. Taxonomic levels are listed in the left margin, with levels which are optional/not applicable to all organisms denoted as such. Rounded boxes shaded in light green correspond to possible taxonomic designations at their respective level, while blue boxes correspond to requirements that an isolate and/or its genome must meet to be assigned that designation. Possible forms which the final taxonomic assignment can take can be found in the gray box at the bottom of the chart.