Bacillus cereus group, also named B. cereus sensu lato (s.l.), contains closely related Gram-positive, facultative aerobe, and endospore-forming bacteria. It includes B. cereus sensu stricto (s.s.), Bacillus thuringiensis, Bacillus anthracis, Bacillus mycoides, Bacillus pseudomycoides, Bacillus weihenstephanensis, and Bacillus cytotoxicus. Recently though, more species have been proposed as new members, such as Bacillus toyonensis (Jimenez et al., 2013), Bacillus wiedmannii (Miller et al., 2016), Bacillus gaemokensis (Jung et al., 2010), Bacillus bingmayongensis (Liu et al., 2014), or Bacillus manliponensis (Jung et al., 2011). The group gathers member species from various ecological niches that display a large spectrum of virulence, from the insect pathogen B. thuringiensis, used worldwide as biopesticide and plant protection agent, to human pathogenic strains of B. anthracis or emetic B. cereus (Agata et al., 1995; Turnbull, 1999). It also spans from psychrotrophic (growth at temperature below 7°C) to thermotolerant (growth up to 53°C) strains (Stenfors Arnesen et al., 2008).

Some strains of B. cereus s.l. have been implicated in food outbreaks associated with emetic or diarrheic syndromes. The former is due to cereulide, a thermostable, protease and pH-resistant dodecadepsipeptide toxin preformed in food, whose genetic determinants are plasmid-borne (Agata et al., 1994; Ehling-Schulz et al., 2004; Hoton et al., 2005). The diarrheal syndrome is presumably caused by one or more potential enterotoxins. In fact, B. cereus can produce several potential enterotoxins, including hemolysin BL (HBL), non-hemolytic enterotoxin (Nhe), enterotoxin FM (EntFM), and cytotoxin K (CytK) (Granum, 1994; Lund et al., 2000; Senesi and Ghelardi, 2010). The latter was originally discovered in a B. cereus strain (NVH 391-98) isolated from an outbreak that caused the death to three elderly persons in France in 1998 (Lund et al., 2000). This strain was further shown to be able to grow up to 53°C. The name B. cytotoxicus was coined for this new thermotolerant “B. cereus” strain (Guinebretière et al., 2013). Due to their ability to thrive at relatively high temperatures, strains of B. cytotoxicus have recently been isolated in geothermal waters (Cavello et al., 2020). Some B. cytotoxicus isolates, including strain NVH 391-98, were shown to over-produce the CytK toxin, as compared to other strains of the same species (Heini et al., 2018). Two variants of cytotoxin K have now been described: CytK-1 and CytK-2. The former is associated with B. cytotoxicus and is more cytotoxic than the second, which is found in some mesophilic strains of B. cereus (Fagerlund et al., 2004).

Bacillus cereus members have been classified into seven phylogenetic groups according to panC typing. All emetic B. cereus and thermotolerant B. cytotoxicus strains fall in groups III and VII, respectively (Guinebretière et al., 2008). These seven groups can also be organized in three genomic clades based on their core and pan-genomes. Interestingly though, the B. cytotoxicus strains form a distinct cluster among the B. cereus members (Fagerlund et al., 2007; Auger et al., 2008; Guinebretière et al., 2008; Bazinet, 2017). This B. cytotoxicus diversity was further supported in recent studies in which isolated strains were classified into four genomic clades (A–D; Stevens et al., 2019). Also, a set of strains isolated from food products were classified into six random amplified polymorphic DNA (RAPD) groups and 11 plasmid profiles (Koné et al., 2019). Although the diversity within the B. cytotoxicus species is now well established, the extent of these genetic and genomic variations remains poorly characterized. The aim of this study was to further explore the genomic and genetic diversity of six B. cytotoxicus strains pertaining to different genetic profiles. Whole genome sequences of these strains were compared to those of other B. cytotoxicus genomes publicly available. Sequences of chromosomal indels, mobile elements, and extrachromosomal molecules were also analyzed in detail.

In a previous study, Koné et al. (2019) classified a collection of 57 B. cytotoxicus strains, isolated from different food products, into six RAPD patterns and 11 plasmid profiles. The six isolates used in the present study were selected on the basis of both their distinct RAPD and/or plasmid profiles. Four strains (E8.1, E17.4, E28.3, and PDT2.12) originated from potato flakes, and two (SM1.1 and SM2.8) were isolated from instant soup. As shown in Table 1, isolates sharing the same RAPD patterns (D for SM1.1 and SM2.8 and A for E8.1 and E28.3) showed different plasmid profiles. As for their plasmid profiles, PDT2.12 had a unique one, while E17.4 shared the same pattern with SM2.8 (Table 1).

As shown in Table 2, the chromosome sizes of the six strains ranged from 4,049,327bp (E17.4) to 4,244,837bp (SM2.8), with GC contents of 35.9–36.0%. As expected from our previous work, all six sequenced isolates contained both small (from none to three) and/or large (from none to two) plasmids with sizes varying from 3,421bp (pE283-3) to 83,570bp (pE81-84). They were all circular, with the exception of pE283-14.

Although it has been demonstrated that B. cytotoxicus strains constituted a remote cluster from the other B. cereus group members (Fagerlund et al., 2007), recent studies have shown that isolates of this species displayed noticeable chromosomal and plasmidial diversities (Koné et al., 2019; Stevens et al., 2019). The first B. cytotoxicus strain was isolated during a severe foodborne outbreak in France that led to three fatalities. Nevertheless, present data on this species suggest a strain-dependent, variable cytotoxicity. Indeed, a recent study suggested that the cytotoxic and fatal potential of this species might be lower than initially thought (Burtscher et al., 2021).

The current study extends insight into this intra-species genetic and genomic diversity. As for the other B. cytotoxicus genomes publicly available (Stevens et al., 2019), the six new strains sequenced in this study have chromosome sizes ranging from 4.1 to 4.2Mb. Also, their genomes contained up to 3.3% of plasmid DNA, with sizes ranging from 3.4 to almost 84kb. Compared to the four recently described clades (i.e., clades A–D, Stevens et al., 2019), two of our strains, E8.1 and E28.3, belong to clade C (Figure 2), that contained the majority of publicly available B. cytotoxicus genomes, while SM1.1 and SM2.8 pertain to clade D, together with strain AFSSA_08CEB44Bac (BioProject: PRJEB14962) isolated in France. The two remaining strains (i.e., PDT2.12 and E17.4) do not belong to previously described clades. This is in line with our recent study showing that these B. cytotoxicus isolates were the sole member of their RAPD patterns, while SM1.1 and SM2.8 displayed the same RAPD pattern (Koné et al., 2019), but differ by 13 chromosomal indels. Of note is the origin of these SM1.1 and SM2.8 strains that were isolated from two Moroccan soups: Harira and Chorba, respectively (Koné et al., 2019). However, it is worth mentioning that these instant soups contained potato-derived products as additives, which might be related to the presence of these bacteria.

Concerning the link between diversity and potential pathogenicity, it has been suggested that strains from clade A (which contains the original and reference type-strain NVH 391-98 and the highly cytotoxic CH_213) and from clade B were more likely to be cytotoxic compared to those of clades C and D (Stevens and Johler, 2020). However, the exact factor(s) and mechanism(s) responsible for the higher cytotoxicity displayed by these B. cytotoxicus strains remain(s) to be further explored.

The functional comparison showed that four out of six sequenced isolates (SM1.1, SM2.8, PDT2.12, and E17.4) possess enzymes implicated in the degradation of inositol, a polyalcohol mainly present in soil and plants. This catabolic operon has been found in environmental dwelling bacteria such as Bacillus subtilis or Klebsiella aerogenes (Yoshida et al., 2008).

As for other B. cereus group members, the B. cytotoxicus reference strain is not able to use galactose (Guinebretière et al., 2013), and analyses of publicly available sequences did not find any galactose degradation enzymatic pathway in their genomes. Surprisingly PDT2.12 contains sequences coding for galactose uptake and degradation enzymes. This hints that this isolate may have thrived in dairy environments and acquired these sequences from other galactose-using bacteria through horizontal gene transfer.

To avoid the taxonomic ambiguity of the B. cereus group, B. cytotoxicus genomes were searched for the presence of virulence genes typically used as markers for the entomopathogenic B. thuringiensis strains, i.e., those encoding crystal proteins or those toxic for insect larvae: the δ-endotoxins cry, cytolytic cyt, and the vegetative insecticidal protein coding genes vip. No entomopathogenic genes were found on the chromosomes or plasmids of the B. cytotoxicus strains.

More than a decade ago, the interest in phages preying on the B. cereus group was revived, bringing forth the questions about their diversity and potential implication in the ecology and adaptability of members of this group (Gillis and Mahillon, 2014a). While some have a lytic life cycle, others are lysogenic and can remain in a dormant state as prophages integrated into the chromosome, into plasmids or as circular/linear elements aka plasmidial prophages (Gillis and Mahillon, 2014a; Piligrimova et al., 2021). Prophages found in B. cereus s.s., B. anthracis, and B. thuringiensis strains have been extensively studied, thanks to the abundance of genomic sequences of these three species, in contrast to B. cytotoxicus genomes, still relatively new to the genomic field. In this study, the genomes of nine B. cytotoxicus strains were analyzed for the presence of prophages, via the online tool PHASTER. Six putative prophages were found to be “intact” on the chromosomes or plasmids. Interestingly, some prophages marked diverse regions on the chromosomes (regions 5, 6, 8, 10, and 11, Figure 1), making them a valuable addition to the bacteria’s genetic pool. For chromosomal and plasmidial prophages, the question of their activity and potential for a lytic cycle is still unanswered and requires further investigation. Nonetheless, for the plasmidial prophages, an added value could be their potential for horizontal genetic transfer, hence their status as “mobilizable.”

A special case is that of the prophage dubbed as B025-like, found as a plasmidial prophage in SM1.1, pSM11-43, and integrated in the chromosome of strain SM2.8. B025 originally isolated from a Listeria strain has the tools required for integration (att sequences and an integrase coding gene) and circularization (cohesive single-strand complementary cos ends). The capacity of a bacteriophage to exist in both integrated and excised forms was previously reported for a Streptococcus pyogenes M1 serotype phage SpyCIM1. The latter forms a chromosomal island integrated at a specific site of the chromosome, while the bacterium is in a stationary growth phase. However, once in exponential growth phase, SpyCIM1 excises from the chromosome and replicates as an independent plasmidial prophage (Nguyen and McShan, 2014; Utter et al., 2014). Another example is the circular plasmidial prophage pLUSID3, recently described in B. thuringiensis strain HER1410, which can also occur as fully integrated within the chromosome causing the disruption of a flagellar key component (Lechuga et al., 2020). Whether these B025-like elements are active or defective prophages and whether some are merely integrated into resident plasmids or are true plasmidial prophages will required further experiments.

Other mobile genetic elements were mined in the B. cytotoxicus genomes, including IS elements. Although no plasmidial IS were recovered, the variability of chromosomal IS elements between the different strains and clades was evident, with NVH 391-98 presenting the lowest number of IS. Surprisingly, 1.62% of the chromosome of SM1.1 is covered by IS, the highest number for the nine strains, which is higher than the average 1.1% of B. thuringiensis, keeping in mind that for the latter species, IS prevalence and diversity are dominantly plasmidial (Fayad et al., 2019).

The plasmidome showed the greatest diversity among the six sequenced B. cytotoxicus strains, including four small plasmids and six large ones, whose prevalence and distribution are summarized in Table 4. While some plasmids were unique to a particular strain, e.g., the 3- and 4-kb plasmids of E28.3, others were more or less similar to plasmids carried by other strains. For the 12-kb plasmids, one set (12a) is identical in SM1.1 and SM2.8, while the other (12b) is identical between SM1.1 and SM2.8 but is closely related to pE174-12 and the 7-kb pBC9801 element of the reference strain NVH 391-98. The main difference between the 7-kb, 12b, and pE174-12 plasmids is the replication protein (Figure 3). However, one of the main shared regions is encoding a Fibronectin type III domain-containing protein. This domain mediates protein–protein interaction and is potentially involved in the correct positioning of a protein’s active domain (Campbell and Spitzfaden, 1994). Although its exact role remains unclear, it is often associated with enzymes involved in the degradation of recalcitrant polysaccharides such as chitin, as shown for the chitinase of B. thuringiensis (Juárez-Hernández et al., 2019).

The most prevalent plasmid, or plasmidial region, is pPDT212-44, to which similarities are found on five other plasmids, in five strains from two clades, and one not placed in a particular clade (Table 4; Figure 4). Finally, the 53-kb conjugative plasmid pE81-53 shows similarities with the other 53-kb elements in previously reported B. cytotoxicus strains, implying their conjugative potential as well (Figure 5; Stevens et al., 2019). In total, the 53-kb element was found in four strains from two different clades, B and C (Table 4).

Another quite interesting feature of the small plasmids is the presence on pSM11-12a of strain SM1.1 (as well as on strain SM2.8) of three genes highly similar to the gakABC locus found in Lactococcus garvieae that code for the three-peptide bacteriocin Garvicin KS (Ovchinnikov et al., 2016). This bacteriocin and homologues found in B. cereus strains are active against several other Gram-positive bacteria including members of the Bacillus, Enterococcus, Listeria, and Streptococcus genera. Intriguingly, the homology extends outside this locus and includes neighboring genes such as putative ABC transporter and recombinase/integrase-like genes, suggesting they are part of a mobile genetic element. The activity of this Garvicin KS-like bacteriocin, also referred to as Cereucin in the case of the B. cereus strains (Ovchinnikov et al., 2016), on other bacteria is currently under investigation.

The aim of this study was to extend the understanding about the intra-species diversity of B. cytotoxicus through WGS and comparative genomic analyses of six isolates. The SNP-based phylo-dendrogram, as well as the ANI, showed that two isolates from instant soup fitted in the clade D, two from potato flakes were classified in clade C, and the last two, also from potato flakes, formed a separated cluster from the other clades. The plasmidial diversity is also in line with previous studies. The presence of sequences coding for inositol degradation in four isolates genomes and galactose uptake and degradation enzymes hinted the ecological niche of B. cytotoxicus. Nevertheless, in order to extend the insight on knowledge about the genetic diversity of B. cytotoxicus and its ecological niche, there is a need to sequence and study more isolates from matrices other than potato products.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/, PRJNA684687.

KK, NF, AG, and JM contributed to conceptualization, validation, methodology, and writing – review and editing. KK, AG, and NF provided software and were involved in investigation. KK, NF, and JM contributed to formal analysis, data curation, writing – original draft preparation, and visualization. JM was involved in resources, supervision, project administration, and funding. All authors contributed to the article and approved the submitted version.

This work was supported by the International Office for Cooperation of the Université catholique de Louvain (UCLouvain; Bursaries to KK and NF), the Research Department of the Communauté française de Belgique (Concerted Research Action, ARC 17/22-084), and the National Fund for Scientific Research (FNRS, Belgium; research grant FNRS-CDR J.0144.20 to JM and research position FNRS 1.B208.16 to AG).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.