# Population structure of Bacillus cereus sensu lato associated with foodborne outbreaks in France between 2004 and 2023

Ksenia Mozhaitseva, Sylvie Pairaud, Olivier Firmesse, Mathilde Bonis

► To cite this version:

Ksenia Mozhaitseva, Sylvie Pairaud, Olivier Firmesse, Mathilde Bonis. Population structure of Bacil- lus cereus sensu lato associated with foodborne outbreaks in France between 2004 and 2023. Food Microbiology, 2025, 133, 10.1016/j.fm.2025.104882. anses-05215966

HAL Id: anses-05215966

https://anses.hal.science/anses-05215966v1

Submitted on 20 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Population structure of *Bacillus cereus sensu lato* associated with foodborne outbreaks in France between 2004 and 2023

Ksenia Mozhaitseva, Sylvie Pairaud, Olivier Firmesse, Mathilde Bonis*

<figure><img src="image_2.png" /></figure>

*Laboratory for Food Safety, French Agency for Food, Environmental and Occupational Health and Safety (ANSES), 94700, Maisons-Alfort, France*

## ARTICLE INFO

**Keywords:**

*Bacillus cereus sensu lato*

*Bacillus cereus sensu stricto*

*Bacillus thuringiensis*

*Bacillus paranthracis*

Foodborne outbreak

Antimicrobial resistance

Bacterial population

## ABSTRACT

*Bacillus cereus sensu lato* (*Bcsl*) is a group of closely related bacterial species known for their resistant spores, enabling them to persist in a dormant state and thereby colonize and adapt across diverse environments. *Bcsl* is known for its harmful impact on human health, producing toxins that cause emetic and diarrheal syndromes or provoking extradigestive infections. Importantly, *Bcsl* is the most frequent confirmed or presumptive causative agent associated with foodborne outbreaks (FBOs) in France. In our study, we assessed the population structure of a large collection of *Bcsl* isolated during FBOs investigation in France between 2004 and 2023, focusing on the association between distinct populations and food categories. Using 294 genomes from 183 FBOs, we applied genomic clustering and phylogenomic analysis and then identified three predominant *Bcsl* populations. *B. cereus sensu stricto* (17.0 %) prevailed in composite dishes, *B. paranthracis* (16.1 %) was positively associated with cereals, and *B. thuringiensis* subsp. *kurstaki* (7.6 %) was predominantly found in vegetable-based salads. Some strains were phylogenetically closely related to clinical isolates, highlighting the need to assess the antibiotic susceptibility of *Bcsl*. Notably, one *Bcsl* clade, *B. cytotoxicus*, lacking beta-lactamase-encoding genes showed a greatly increased sensitivity to ampicillin than other *Bcsl* considered to be naturally resistant to beta-lactams. Additionally, various strains from distinct populations showed reduced susceptibility to macrolides and cyclins. Finally, accurately differentiated populations will be used in further epidemiological studies and in dose-response modeling.

## 1. Introduction

*Bacillus cereus sensu lato* (*Bcsl*), or the *Bacillus cereus* group, comprises numerous validly and effectively published species (Parte et al., 2020). These bacteria are Gram-positive, aerobic or facultatively anaerobic, spore-forming, and primarily telluric and rhizospheric (Vilain et al., 2006). The spores of *Bcsl* are resistant to heat, UV radiation, acid, and desiccation, enabling them to persist in a dormant state and thereby contaminate diverse environments (Ehling-Schulz et al., 2019). In addition to soil, *Bcsl* has been isolated from plants, seawater, food products, animals, and humans (Schnepf et al., 1998; Guinebretière et al., 2013; Miller et al., 2016; Liu et al., 2017; Glasset et al., 2018; Tohya et al., 2021; Muigg et al., 2022; Makuwa et al., 2023; Robas Mora et al., 2023). *Bcsl* includes species known for their harmful impact on human health.

It is the most frequent confirmed or presumptive causative agent associated with foodborne outbreaks (FBOs) in France (Santé Publique France. Bulletin, 2024), and the second most frequently reported

causative agent in Europe (European Food Safety Authority, European Centre for Disease Prevention and Control, 2024). *Bcsl* can produce various toxins, resulting in emetic or diarrheal syndromes. The emetic syndrome is caused by cereulide, a heat-stable, non-ribosomal cyclic peptide, whose synthesis is directed by the *ces* operon located in the pCER270 megaplasmid. In the most severe cases, ingestion of this toxin, produced in food, can lead to multiple organ failure and death (Rouzeau-Szynalski et al., 2020). Cereulide-producing strains are commonly referred to as emetic *B. cereus* and are predominantly found in specific *B. paranthracis* clades (Carroll et al., 2020b; Carroll et al., 2022b; Frentzel et al., 2024). However, they are also occasionally identified in other *Bcsl* (Mei et al., 2014; Pheepakpraw et al., 2023; Frentzel et al., 2024). The diarrheal syndrome is associated with the pore-forming enterotoxins: hemolysin BL (Hbl), non-hemolytic enterotoxin (Nhe), and cytotoxin K (CytK). Hbl and Nhe are tripartite toxins encoded by the *hblCDA* and *nheABC* operons, respectively (Tuipulotu et al., 2021). Cytotoxin K is encoded by the *cytK* gene existing in two variants: *cytK-1* and *cytK-2*. CytK-1 shows higher cytotoxicity and is

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: mathilde.bonis@anses.fr (M. Bonis).</footnote>
restricted to *B. cytotoxicus* (Guinebretière et al., 2013), while CytK-2 can be produced by other *Bcsl* (Fagerlund et al., 2004). These enterotoxins, produced in the digestive tract, damage intestinal epithelial cells by forming pores leading to osmotic lysis and, ultimately, diarrhea. Although less characterized, other virulence factors, such as sphingomyelinase, enterotoxin FM, hemolysins I, II and III, metalloproteases, and phospholipases complement the action of the *Bcsl* enterotoxins (Tuipulotu et al., 2021).

The members of the *B. cereus* group are also able to provoke extra-digestive infections, such as anthrax (*B. anthracis*). In addition, *Bcsl* has been shown to disseminate in hospital environments contaminating medical devices, equipment, and linen, and cause various local and systemic infections, particularly in immuno-compromised individuals (Glasset et al., 2018). Certain species, such as B. sanguinis and B. par-amobilis, were originally isolated from the blood of hospitalized patients (Tohya et al., 2021), and B. basilensis was first recovered from a wound swab under clinical conditions (Muigg et al., 2022). Moreover, in B. mobilis, originally isolated from ocean sediment (Liu et al., 2017), one strain was identified as clinically relevant and linked to polymicrobial catheter infection and enteritis (Muigg et al., 2022). This reinforces the need to properly define the antibiotic susceptibility pattern of *Bcsl*. To date, *Bcsl* has been described as naturally resistant to beta-lactams, while resistance to other families of antibiotics were described more sporadically (Glasset et al., 2018; Bianco et al., 2021; Jung et al., 2022; Rajalingam et al., 2022; Mohammadi et al., 2023).

Various *Bcsl* are known for their commercial use, in particular, *Bacillus thuringiensis* (*Bt*) which can produce delta-endotoxins (Cry, Cyt, and Vip) with insecticidal activity (Schnepf et al., 1998). Several *Bt* strains, primarily those belonging to subspecies *aizawai* (*Bta*) and *kurstaki* (*Btk*), have been commercialized and are widely used for pest control (Jallouli et al., 2020). Besides, certain strains of *B. toyonensis* are used as probiotics in animal feed (Pandith and Kumar, 2024). An isolate of another species, *B. mycoides*, is used as a fungicide for plant protection (Jacobsen et al., 2007).

Since *Bcsl* is of particular importance for FBO surveillance, it is crucial to accurately differentiate *Bcsl* strains and identify biologically meaningful groups to further associate them with relevant data, such as their prevalence across various food categories.

However, the *B. cereus* group is known for the challenges in taxonomy and proper definition of population boundaries due to its high genome plasticity (Ehling-Schulz et al., 2010; Böhm et al., 2015; Bazinet, 2017) and elevated genetic similarity between some species (Liu et al., 2015; Bazinet, 2017; Torres Manno et al., 2020). Historically, *Bcsl* members were classified into different species primarily based on phenotypic traits, such as colony morphology (Flügge, 1886), growth temperature specificity (Lechner et al., 1998), or association with symptoms and diseases (Mortimer and McCann, 1974). Later, other criteria included presence of virulence plasmids, similarity measured using DNA-DNA hybridization techniques, phylogenetic groups and sequence types based on comparison of single locus and multiple loci (SLST and MLST), respectively (Carroll et al., 2022a). When applying the SLST approach, partial sequence of the gene encoding pantoate-β-alanine ligase (*panC*) allows a relatively high resolution to differentiate the major phylogenetic groups of *Bcsl*, I to VII (Guinebretière et al., 2008). To date, *panC* phylogenetic groups I, V, and VII are represented by *B. pseudomycoides*, *B. toyonensis* and *B. cytotoxicus*, respectively. Group IV includes *B. cereus sensu stricto* (*Bcss*) and *Bt*. Groups II, III, and VI are the most diverse and comprise several species (Carroll et al., 2020a). Among MLST schemes, the one based on seven housekeeping genes (Priest et al., 2004) is the most widely used. For example, using this scheme, certain sequence types (STs) were found to be frequently associated with foodborne outbreaks, such as emetic strains of *B. paranahracis* belonging to ST26 (Frentzel et al., 2024). The study of Tong et al. (2024) found that within *Bcss*, ST24 and ST142, both associated with foodborne illness, were the most prevalent STs in Europe.

Finer characterization of *Bcsl* is feasible thanks to whole genome sequencing data (Liu et al., 2015; Bazinet, 2017; Torres Manno et al., 2020). Currently, the assignment of *Bcsl* strains to genospecies is based on the calculation of the Average Nucleotide Identity (ANI) index, the percentage of nucleotide identity in the alignment between query and reference (e.g. type strain) genomes. Commonly, an ANI equal to or higher than 95 % reflects natural species boundaries (Jain et al., 2018). However, this 95 % threshold is inadequate for the *B. cereus* group, leading to ambiguous genospecies assignment results, as discussed in Carroll et al. (2020b, 2021). To avoid overlapping genospecies within *Bcsl*, the authors proposed an ANI threshold of 92.5 %. Using this threshold, both *Bcss* and *Bt* are grouped into the same *Bcss* genospecies, while all species belonging to *panC* groups II (except for *B. luti*) and III are grouped into the *B. mosaicus* genospecies. Any strain carrying a plasmid that harbors Cry/Cyt/Vip toxin, emetic toxin and/or anthrax toxin encoding genes is labeled with the corresponding biovar name (Thuringiensis, Emeticus, and/or Anthracis, respectively), appended to the genospecies name. This nomenclature is useful for clinicians as it helps avoid ambiguous genospecies assignment and enables the rapid identification of genotypes with clinical and industrial relevance. However, it is less suited for finer resolution, particularly for *panC* groups III and IV. *Bt*, widely used as a biopesticide, can be confused with the phylogenetically closely related *Bcss*, a known food poisoning agent. The latter may carry the plasmids harboring insecticidal toxin genes acquired via horizontal transfer (Battisti et al., 1985). Furthermore, these two species show similar phenotypic characteristics (van Tongeren et al., 2014), highly similar chromosomal structure (Carlson et al., 1996), and an ANI greater than 96 % (Torres Manno et al., 2020). Therefore, it is of particular importance to accurately distinguish *Bt* from *Bcss*, given increasing data on the *Bt* presence in food products, its association with foodborne outbreaks (Frentzel et al., 2020; Bonis et al., 2021; Biggel et al., 2022a, 2022b) and nosocomial infections (Kuroki et al., 2009; Bianco et al., 2021; Barakat et al., 2024).

*Bcsl* is characterized by a dynamic pangenome due to its high genome plasticity (Ehling-Schulz et al., 2010; Böhm et al., 2015; Bazinet, 2017). The pangenome of *Bcsl* includes a core genome, which contains genes shared by 95-100 % of strains within a bacterial species, and accessory genome, which consists of genes found in only some strains. Using 146 *Bcsl* genome assemblies, Bazinet (2017) estimated its pangenome at ~60,000 genes, while the core genome at ~600 genes, thereby suggesting high adaptability of *Bcsl* resulting from the acquisition of accessory genes. *Bcsl* accessory genes are enriched in various functions including those linked to antimicrobial resistance and virulence (White et al., 2022). Furthermore, Böhm et al. (2015) demonstrated that toxin-encoding accessory genes, such as *hbl* and *cytK*, were horizontally transferred between *Bcsl* clades. The authors concluded that this could enable non-pathogenic strains to acquire virulence factors. This might have consequences for risk assessment of *Bcsl*. As a result, an accurate characterization of a particular *Bcsl* collection of interest requires multiple complementary methods, including reverse ecology approach, phylogenomics, MLST, and accessory gene content analysis. The aim of this study was to resolve the population structure of *Bcsl* isolated from foodborne outbreaks in France between 2004 and 2023, with a focus on the relationship between predominant populations and specific food categories. We examined a large *Bcsl* collection, consisting of 294 *Bcsl* genomes from 183 FBOs. To avoid the challenges in *Bcsl* taxonomy and accurately define species and population boundaries, we first addressed the reverse ecology approach and applied an appropriate tool (Arevalo et al., 2019). This helped us differentiate the samples into biologically meaningful units and perform further characterization of populations of interest by performing phylogenomic, MLST, admixture and accessory gene content analyses. We identified three predominant FBO-associated *Bcsl* populations in France, each statistically associated with a specific food category. We consider finely differentiated *Bcsl* populations to be a key variable to perform further epidemiological analyses and dose-response modeling. In this perspective, these findings
may be considered to improve the surveillance of *Bcsl*.

## 2. Materials and methods

### 2.1. Genome assemblies

We analyzed a collection of 294 draft genomes of *Bcsl* (Supplementary Table S1) sampled during foodborne outbreaks in France between 2004 and 2023 and stored in the Food Safety Laboratory strain bank of the French Agency for Food, Environmental and Occupational Health and Safety (ANSES). The corresponding raw reads are available at the NCBI platform under the projects PRJNA1209675, PRJNA781790, and PRJNA547495. The DNA extraction technique and the sequencing protocol are described in Bonis et al. (2021). The genomes were assembled using the in-house *bacflow* pipeline, as described in Fichant et al. (2024). The assembly size varied from 4.1 to 7.0 Mbp, the N50 statistic ranged from 15.9 to 936 kbp, and the L50 values ranged from 2 to 91. The assembled genomes of 23 clinical and blood culture isolates were retrieved from the MicroBIG-E browser (https://www.ncbi.nlm.nih.gov/pathogens/microbigge/, accessed in June 2024). The data were filtered by location (France) and isolation source (blood, patient, intravascular catheter) (Supplementary Table S2).

### 2.2. Species tree

An unrooted phylogenetic tree was constructed based on 1723 protein sequences of single-copy orthologous genes identified across the *Bcsl* genomes using OrthoFinder v.2.5.5 (Emms and Kelly, 2015). Multiple protein sequence alignment was performed by MUSCLE v.August 3, 1551 (Edgar, 2004). The resulting 'SpeciesTreeAlignment.fa' file was used to infer a maximum likelihood phylogeny in IQ-TREE v.2.2.6 (Minh et al., 2020). The best protein substitution model identified by ModelFinder (Kalyaanamoorthy et al., 2017) (-m MFP tag in IQ-TREE) was JTT + F + I + R10 (Jones et al., 1992; Yang, 1995; Soubrier et al., 2012). The number of replicates for ultrafast bootstrap approximation (Minh et al., 2013; Hoang et al., 2018) for the branch support assessment was set to 1000 (-B 1000 tag in IQ-TREE).

### 2.3. Taxonomic assignment, panC group, and Bt typing

Taxonomic assignment of *Bcsl* genomes was performed using the framework of Carroll et al. (2020b) implemented in Btyper3 v.3.4.0 (Carroll et al., 2020a). For each genome, the tool identified the genospecies, the closest type strain, and the PubMLST (https://pubmlst.org) sequence type. The samples were further assigned to one of the *panC* (I to VII) phylogenetic groups proposed by Guinebretière et al. (2013) using a publicly available tool (https://toolcereusid.shinyapps.io/Bcereus/; accessed in March 2024). In addition, *Bta* and *Btk* were identified based on the presence of *Bt*-specific genes following the workflow described in Fichant et al. (2022) and implemented in *Bt typing* (https://github.com/afelten-Anses/Bt Typing; accessed in March 2024).

### 2.4. Genome clustering

The clustering of *Bcsl* genomes was performed using PopCOGenT (Arevalo et al., 2019) (https://github.com/philarevalo/PopCOGenT, accessed February 2024), designed to identify bacterial populations by calculating a distribution of identical DNA sequence lengths between pairs of genomes. PopCOGenT first reconstructs networks (clusters) of similar genomes and further resolves subclusters representing putative populations.

### 2.5. Pangenome examination

The *Bcsl* genomes were annotated with Bakta v.1.9.3 (Schwengers

et al., 2021) using the Bakta database v.5.1 (January 2024). The resulting .gff3 annotations were processed with Panaroo v.1.4.1 (Tonkin-Hill et al., 2020) with default parameters. The pangenome was reconstructed at two levels: (i) across all members of *Bcsl* and (ii) within a specific PopCOGenT cluster of interest. Further hierarchical clustering of the *Bcsl* genomes was based on gene presence-absence and performed using the gene presence absence Rtab output matrix from Panaro. Pairwise Jaccard distances were calculated by applying the dist (method = "binary") R function. The resulting distance matrix was used in the pvclust (method.hclust = "complete", nboot = 1000) package v.2.2.0 (Suzuki et al., 2019) implemented in R v.4.3.2 (R Core Team, 2023). The genomes sharing similar accessory gene content were clustered into the Groups. The Group was considered significant at approximately unbiased p-value >95 % and bootstrap probability >95 %. Thereafter, candidate Group-specific genes were identified by Scoary v.1.6.16 (Brynildsrud et al., 2016), applying the same parameters as Fichant et al. (2022): with thresholds of 80 % sensitivity, 80 % specificity, and a p-value cutoff of 0.05 corrected for multiple testing by Benjamini-Hochberg method. Further Cluster of Orthologous Groups (COG) enrichment analysis was performed in MicrobiomeProfiler (Wu et al., 2021) at two levels: (i) PopCOGenT cluster against the entire background *Bcsl* COG set, and (ii) Groups found within the PopCOGenT cluster of interest against the entire background cluster COG set. Significant COGs were identified at an adjusted p-value cutoff set to 0.05 (Benjamini-Hochberg correction).

### 2.6. Population structure

We carried out a population structure analysis for each composite PopCOGenT cluster consisting of several subclusters. First, we applied fastGEAR (Mostowy et al., 2017) which uses hidden Markov model to assess shared ancestry and reveal putative recombining regions within homologous regions of bacterial genomes. The tool defines lineages as groups of samples that share identical sequences at least 50 % of the sites. Recombination events were inferred by identifying regions with a log Bayes factor (logBF) > 0 indicating a higher probability of recombination than mutation. We focused on both recent and ancestral recombination. Ancestral recombination was defined by fastGEAR as putative recombination event affecting all samples in the lineage, while recent recombination referred to events occurring in some but not all isolates. As input, we used a filtered alignment of a set of core genes reconstructed by Panaroo v.1.4.1 (Tonkin-Hill et al., 2020) for each PopGOGenT cluster. Ambiguous sites and indels were excluded from the alignment using snp-sites v.2.5.1 (Page et al., 2016). Second, we validated results by assessing recombination with ClonalFrameML (Didelot and Wilson, 2015) which employs a maximum likelihood method. The program estimates an average frequency of recombination over mutation (R/θ), an average number of substitutions (δv) per recombination event, and an overall impact of recombination on phylogeny (r/m = R/θ × δv) in a population of interest. A guide phylogenetic tree was reconstructed for each cluster with IQ-TREE v.2.2.6 (Minh et al., 2020) using the filtered alignment of core genes, the best nucleotide substitution model identified by ModelFinder (Kalyaanamoorthy et al., 2017), and 1000 replicates for ultrafast bootstrap approximation ((Minh et al., 2013; Hoang et al., 2018). Thereafter, the tree topology was refined by excluding regions affected by recombination.

### 2.7. Identification of toxin, antimicrobial and heavy metal resistance genes

The cesABCD, nheABC and hblCDA operons, cytK-1 and cytK-2 genes, and the genes encoding insecticidal toxins (Cry, Cyt, and Vip) were detected using Btyper3 v. 3.4.0 (Carroll et al., 2020a) with default parameters. BTyper3 also automatically searched for the *B. anthracis* virulence factors (*pagA*, *lef*, *cya*). The *hlyII* gene was directly searched in the genome annotations. Antimicrobial and heavy metal resistance
genes were identified by AMRFinderPlus v.3.12.8 applying default parameters and using the reference protein database v.2024-01-31.1 (Feldgarden et al., 2021).

## 2.8. Antibiotic susceptibility testing

Antibiotic susceptibility of the *Bcsl* isolates (Supplementary Tables S1 and S4) was assessed using the Kirby-Bauer disk diffusion method in accordance with the European Committee on Antimicrobial Susceptibility Testing (EUCAST) methodology (method v13.0 and breakpoint tables for interpretation v15.0; https://www.eucast.org; accessed in January 2025). Briefly, after isolation on trypticase soy agar-yeast extract (TSA-YE, bioMérieux) medium for 18h at 30 °C, bacteria were resuspended in peptone salt broth (bioMérieux) at 0.5 McFarland turbidity and inoculated onto Mueller Hinton agar (Condalab) plates. Antibiotic disks (Oxoid) containing ampicillin (10 µg, AMP10), erythromycin (15 µg, E15), tetracyclin (30 µg, TE30) and vancomycin (5 µg, VA5) were then applied onto inoculates. Plates were incubated for 18h at 35 °C before measuring the inhibition diameters.

## 2.9. Proportions of FBO-related *Bcsl* populations

To minimize the risk that the set of 294 available sequenced samples might not accurately reflect the true proportions of the FBO-associated *Bcsl* populations, we filtered the dataset as follows. First, we reviewed all ANSES FBO reports (n = 560) indicating the presence of *Bcsl*, including both sequenced and non-sequenced samples. The reports provided information on (i) the *panC* phylogenetic group determined by PCR, and (ii) the microscopic observation of insecticidal toxin crystals (*Bt* phenotype). Second, based on the reports, we calculated the proportion of each *panC* group: I – 0 %, II – 15 %, III – 38 %, IV – 35 %, V – 5 %, VI – 5 %, and VII – 2 %. Within *panC* group IV, we identified the *Bt* phenotype based on its ability to produce insecticidal crystals. The *Bcsl*: *Bt* ratio was estimated at 2:1. Third, we applied the bootstrap technique by generating 100 random sample sets in R v.4.3.2 (R Core Team, 2023) with respect to the obtained proportions, allowing no more than 1 % of deviation. In each random sample set, the maximum possible number of samples, while respecting necessary proportions, was 147. Finally, from the generated sets, we retrieved the average percentage for each *Bcsl* population.

## 2.10. Association of *Bcsl* populations with food categories

After estimating proportions of the *Bcsl* populations, we then examined their distribution across various food categories. Information on contaminated products (Supplementary Table S1) was available in the FBO reports. Food items were classified using both the FoodEx2 (exposure hierarchy: Level 1 and 2) (European Food Safety Authority, 2015) and Food and Agriculture Organization (FAO) (Level 1 and 2) (Leclercq et al., 2019) classifications. Mean percentages of the *Bcsl* populations were calculated for each food category, based on the previously simulated datasets (see above). To determine whether a specific *Bcsl* population was predominantly associated with a particular food category, we performed a Fisher's exact test (fisher.test function) using R v.4.3.2 (R Core Team, 2023). For each food category, we constructed a contingency table comparing mean counts of the *Bcsl* population of interest with the mean counts of the entire *Bcsl* group (background). Associations were considered significant at *p* < 0.05.

## 2.11. Data visualization

Phylogenetic trees and plots were designed using ggtree v.3.5.1 (Yu et al., 2017; Xu et al., 2021; Yu, 2022), ape v.5.8 (Parte et al., 2020) and ggplot2 v.3.10.1 (Wickham, 2016) packages implemented in R v.4.3.2 (R Core Team, 2023).

# 3. Results

## 3.1. *Bcsl* phylogeny and clustering

The studied *Bcsl* isolates were assigned to seven *panC* phylogenetic groups and divided into 14 PopCOGenT (0–13) clusters and 12 unclustered samples (Fig. 1, Table 1, Supplementary Table S1). A genome was considered unclustered if it represented a single sample within the cluster, assigned to a specific *Bcsl* species. For eleven clusters, PopCOGenT did not suggest any division into subclusters. We identified cluster 2 (n = 15, *B. cytotoxicus*), cluster 3 (n = 5, *B. mobilis*), cluster 6 (n = 11, *B. mycoides*), cluster 7 (n = 3, *B. luti*), cluster 8 (n = 9, *B. toyonensis*), cluster 12 (n = 2, *B. paramobilis*), and cluster 13 (n = 2, *B. wiedmannii*). Both cluster 4 (n = 11) and cluster 11 (n = 2) showed the highest ANI with the *B. anthracis* type strain, with respective ANI 97.2 to 98.7 and 94.3. The samples from both clusters did not possess the *pagA*, *lef*, and *cya* anthrax-associated genes. Two other clusters, 9 (n = 8, *Bt*) and 10 (n = 2, *Bt*), showed the highest ANI with the *Bt* Berliner type strain: 96.0–96.4 and 97.3–98.3 ANI, respectively. By contrast, three PopCOGenT clusters demonstrated composite structure and contained several subclusters (populations). Cluster 0 (n = 75) consisted of nine PopCOGenT subclusters (0.0–0.8) containing samples exhibited the highest ANI with either the *B. paranthracis*, *B. pacificus* or *B. tropicus* type strains. Therefore, the cluster was named the *B. paranthracis* + *B. pacificus* + *B. tropicus* cluster. Another composite cluster, cluster 1 (n = 113), comprised samples closely related to the *Bcsl* Frankland and Frankland type strain and divided into subclusters 1.0 and 1.1. The cluster was named *Bcsl* + *Bt* as additional analysis revealed that one of the subclusters included several *Bt* strains (see below). Finally, cluster 5 (n = 24) comprised samples closely related to either the *B. wiedmanni* or “*B. pretiosus*” type strains and therefore named the *B. wiedmanni* + “*B. pretiosus*” cluster. “*B. pretiosus*” is placed in quotation marks because it is an effectively published name (Robas Mora et al., 2023) but not validly published under the rules of bacterial nomenclature (Oren et al., 2023). The cluster was divided into five subclusters, 5.0 to 5.4. Subsequently, finer population structure analysis was performed on these three composite clusters.

## 3.2. Population structure of *B. paranthracis* + *B. pacificus* + *B. tropicus* cluster 0

The first composite PopCOGenT cluster, cluster 0, consisted of nine distinct subclusters, 0.0 to 0.8 (Fig. 2A, Table 1). Six of these subclusters (0.0, 0.3, 0.4, 0.5, 0.6, and 0.7) were closely related in terms of ANI to the *B. paranthracis* type strain. Among them, subcluster 0.3 was the largest and included ST26, the most dominant sequence type. Subclusters 0.1 and 0.2 represented two distinct *B. pacificus* populations, while subcluster 0.8 corresponded to a *B. tropicus* population. Additional phylogenetic analysis revealed that the various clinical *B. paranthracis* samples did not form a distinct population. These clinical samples were phylogenetically close to either subclusters 0.3, 0.4, or 0.5. The most prevalent ST within the clinical samples was ST26 (11 out of 16 samples). Subsequently, fastGEAR detected lineages and admixture blocks indicating either shared ancestry or putative recombination (Fig. 2B and C, Table 1). The first fastGEAR lineage grouped the PopCOGenT subclusters 0.0, 0.3, 0.4 and 0.5. The second lineage consisted of subclusters 0.6 and 0.7. The third lineage consisted of subclusters 0.1 and 0.2. The fourth lineage corresponded to subcluster 0.8. Next, fastGEAR assessed putative recent recombination events. The highest total length of recombining regions was detected in the sample 18SBCL174 (subcluster 0.7, 1.11 Mbp), followed by 11CEB28 (0.8; 936.1 kbp), 22SBCL0891_1 (0.7; 831.6 kbp), and 19SBCL1079 (0.4; 456.6 kbp) (Fig. 2C). However, further analysis with ClonalFrameML showed the sample 11CEB28 (subcluster 0.8) to exhibit the highest sum length affected by recombination, estimated at 91.5 kbp, followed by 14SBCL19 (0.8; 30.2 kbp), 22SBCL1010_3 (0.6; 26.9 kbp), 22SBCL0891_1 (0.7; 25.8 kbp), and
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Maximum likelihood phylogenetic tree of *Bcsl* reconstructed using single-copy orthologous protein sequences. Branch lengths are presented in units of the expected number of amino acid substitutions per site. Bootstrap values are shown on internal nodes. Each sample label is colored according to the PopCOGenT cluster and corresponds to networks designed in Gephi v. 0.10.1. Each node corresponds to one sample or clonal complex, i.e. genomes showing high similarity (Arevalo et al., 2019). **A:** *panC* phylogenetic groups. **B:** *genomospecies* identified using the framework of Carroll et al. (2020a, 2020b). **C:** Toxin-encoding genes and operons: emetic toxin (*cesABCD*), insecticidal toxins (*cry/cyt/vip*), cytotoxin K-1 (*cytK-1*), cytotoxin K-2 (*cytK-2*), hemolysin BL (*hblCDA*), hemolysin II (*hlyII*), non-hemolytic enterotoxin (*nheABC*). **D:** Antimicrobial and heavy metal resistance genes: arsenic resistance operon (*ars*), beta-lactamase (*bla1*, *bla2*, *blaIII*), mercury resistance operon (*mer*), erythromycin resistance (*mphL*, *mphM*), tetracycline resistance (*tet(45)*, *tet(L)*), regulatory part of the *van* operon (*vanRS*).</figcaption></figure>

18SBCL174 (0.7; 24.8 kbp) (Fig. 2A). For the entire cluster, the R/θ ratio was estimated at 0.75, indicating recombination to occur about 1.3 times less than mutation. However, the average number of substitutions introduced per recombination event (δv) was 7.7. Thus, the r/m parameter was 5.8. This indicates that, in this cluster, the recombination process overall led to 5.8 times more substitutions than the mutation process.

Gene content clustering revealed that the samples from *B. pacificus* subcluster 0.2 shared a similar set of accessory genes (Group 70) (Fig. 2A, Supplementary Fig. S1). Likewise, the analysis revealed a similar accessory genome within *B. paranthracis* subcluster 0.4 (Group 72). Notably, Group 76, which included genomes from *B. paranthracis* subclusters 0.0 and 0.5, was enriched in accessory genes associated with the ‘Cell wall/membrane/envelope biogenesis’ function (M COG) (Supplementary Table S3). Additionally, within *B. paranthracis* ST26 (subcluster 0.3), we identified two distinct groups, Group 33 and Group 41, each characterized by its own unique set of accessory genes.

### 3.3. Population structure of *Bcss* + *Bt* cluster 1

The second composite PopCOGenT cluster, which included the highest number of samples (n = 113), composed of two subclusters: *Bcss* subcluster 1.0 and *Bt* subcluster 1.1 (Fig. 3A, Table 1). An additional application of the workflow to identify *Bt*-specific markers (Fichant et al., 2022) revealed a finer structure within the *Bt* subcluster. This analysis identified *Bt* subsp. *aizawai* (*Bta*, ST15), *Bt* subsp. *kurstaki* (*Btk*, ST8), two *Bt* ST13 samples (1.1 *Bt*), and four genomes in which no *Bt*-specific markers described in Fichant et al. (2022) were detected (designated hereafter as 1.1 “non-*Bt*”). Further fastGEAR analysis revealed admixture blocks shared between the 1.1 “non-*Bt*” genomes and *Bta*, *Btk*, or 1.1 *Bt* ST13 (Fig. 3B). ClonalFrameML similarly detected

recombination in the 1.1 “non-*Bt*” samples and in the *Bcss* subcluster 1.0. The maximum length of recombining regions was detected for the sample 15SBCL85 (1.1 “non-*Bt*”; 46.8 kbp), followed by 08CEB85 (*Bcss*; 21.9 kbp), 19SBCL1023 (1.1 “non-*Bt*”; 14.5 kbp) and 17SBCL1086 (*Bcss*; 14.5 kbp). For the entire *Bcss* + *Bt* cluster 1, the average R/θ, δv and r/m parameters were estimated at 0.59, 7.3 and 4.3, respectively. Additionally, the reconstructed phylogeny revealed that seven French clinical *Bcss* samples did not form a distinct population but were phylogenetically closely related to various FBO-related strains (Fig. 3A).

Gene content analysis revealed several groups of genomes with similar accessory genes. Within the *Bt* subcluster 1.1, we identified the *Bta* (Group 85) and *Btk* (Group 92) significant groups (Fig. 3A, Supplementary Fig. S2). Within the *Bcss* subcluster, the largest significant group (Group 95) comprised 11 samples corresponding to *Bcss* ST24. Interestingly, the 1.1 “non-*Bt*” samples showed greater similarity to *Bcss* than to *Bta* or *Btk* in terms of gene content. Furthermore, enrichment analysis performed for the significant Groups revealed that *Bta* and *Bcss* ST24 each possessed a specific set of accessory genes enriched in the ‘Mobilome: prophages, transposons’ functional category (X COG) (Supplementary Table S3).

### 3.4. Population structure of *B. Wiedemanni* + “*B. pretiosus*” cluster 5

The third analyzed composite cluster, cluster 5, comprised two species, *B. wiedemanni* and “*B. pretiosus*”. Within this cluster, PopCOGenT identified five distinct subclusters (populations), 5.0 to 5.4 (Fig. 4A, Table 1). Admixture analysis, performed by fastGEAR, identified three distinct lineages (Fig. 4B–Table 1). The first lineage grouped subclusters 5.0 and 5.3, the second grouped subclusters 5.1 and 5.2, and the third corresponded to subcluster 5.4. In the core genome alignment (Fig. 4B), we identified several admixture blocks shared between lineages (5.0 +
<table><caption>**Table 1**
Clustering of *Bacillus cereus sensu lato* associated with foodborne outbreaks in France.</caption><thead><tr><th>PopCOGenT<br>cluster number</th><th>Number<br>of<br>samples<br>per<br>cluster</th><th>Cluster name<br>(closest type strain)</th><th>PopCOGenT subcluster<br>(population) number<br>and name</th><th>Number of<br>samples<br>per<br>subcluster</th><th>fastGEAR<br>lineage</th><th>Bt typing<br>(Fichant et al.,<br>2022)</th></tr></thead><tbody><tr><td rowspan="10">0</td><td rowspan="10">75</td><td rowspan="10">B. paranthracis +<br>B. pacificus +<br>B. tropicus</td><td>0.0 (B. paranthracis)</td><td>5</td><td rowspan="3">(0.0 + 0.3 +<br>0.4 + 0.5)</td><td>-</td></tr><tr><td>0.3 (B. paranthracis)</td><td>36</td><td>-</td></tr><tr><td>0.4 (B. paranthracis)</td><td>10</td><td>-</td></tr><tr><td>0.5 (B. paranthracis)</td><td>3</td><td rowspan="2">(0.6 + 0.7)</td><td>-</td></tr><tr><td>0.6 (B. paranthracis)</td><td>4</td><td>-</td></tr><tr><td>0.7(B. paranthracis)</td><td>2</td><td rowspan="2">(0.1 + 0.2)</td><td>-</td></tr><tr><td>0.1 (B. pacificus)</td><td>4</td><td>-</td></tr><tr><td>0.2 (B. pacificus)</td><td>8</td><td rowspan="2">0.8</td><td>-</td></tr><tr><td>0.8 (B. tropicus)</td><td>3</td><td>-</td></tr><tr><td rowspan="10">1</td><td rowspan="10">113</td><td rowspan="10">Bcss (B. cereus sensu<br>stricto) + Bt (B.<br>thuringiensis)</td><td>1.0 (Bc)</td><td>28</td><td>1.0</td><td>-</td></tr><tr><td rowspan="4">1.1 (Bt)</td><td>26</td><td>1.1 Bta *</td><td>Bta</td></tr><tr><td>53</td><td>1.1 Btk</td><td>Btk</td></tr><tr><td>2</td><td>1.1 Bt<br>(ST13)</td><td>Other Bt<br>lineage</td></tr><tr><td>4</td><td>1.1<br>« non-Bt »</td><td>No Bt markers<br>found</td></tr><tr><td>2</td><td>15</td><td>B. cytotoxicus</td><td>-</td><td>-</td><td>-</td></tr><tr><td>3</td><td>5</td><td>B. mobilis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>4</td><td>11</td><td>B. anthracis</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="10">5</td><td rowspan="10">24</td><td rowspan="10">B. wiedmannii +<br>"B. pretiosus"</td><td>5.0</td><td>2</td><td rowspan="2">(5.0 + 5.3)</td><td>-</td></tr><tr><td>5.3</td><td>5</td><td>-</td></tr><tr><td>5.1</td><td>11</td><td rowspan="2">(5.1 + 5.2)</td><td>-</td></tr><tr><td>5.2</td><td>2</td><td>-</td></tr><tr><td>5.4</td><td>4</td><td>5.4</td><td>-</td></tr><tr><td>6</td><td>11</td><td>B. mycoides</td><td>-</td><td>-</td><td>-</td></tr><tr><td>7</td><td>3</td><td>B. luti</td><td>-</td><td>-</td><td>-</td></tr><tr><td>8</td><td>9</td><td>B. toyonensis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>9</td><td>8</td><td>Bt</td><td>-</td><td>-</td><td>No Bt markers<br>found</td></tr><tr><td>10</td><td>2</td><td>Bt</td><td>-</td><td>-</td><td>No Bt markers<br>found</td></tr><tr><td>11</td><td>2</td><td>B. anthracis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>12</td><td>2</td><td>B. paramobilis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>13</td><td>2</td><td>B. wiedmannii</td><td>-</td><td>-</td><td>-</td></tr><tr><td>04CEB13 **</td><td>1</td><td>B. hominis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>09CEB32</td><td>1</td><td>B. mobilis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>09CEB56</td><td>1</td><td>B. toyonensis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>10CEB25</td><td>1</td><td>B. mobilis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>10CEB26</td><td>1</td><td>B. mobilis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>10CEB29</td><td>1</td><td>B. albus</td><td>-</td><td>-</td><td>-</td></tr><tr><td>14SBCL714</td><td>1</td><td>B. wiedmannii</td><td>-</td><td>-</td><td>-</td></tr><tr><td>15SBCL1065</td><td>1</td><td>B. pseudomycoides</td><td>-</td><td>-</td><td>-</td></tr><tr><td>16SBCL1267</td><td>1</td><td>B. mobilis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>20SBCL79</td><td>1</td><td>B.<br>weihenstephanensis</td><td>-</td><td>-</td><td>-</td></tr><tr><td>22SBCL0394_1</td><td>1</td><td>B. proteolyticus</td><td>-</td><td>-</td><td>-</td></tr><tr><td>22SBCL0404_5</td><td>1</td><td>B. mobilis</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<footnote>* When performing fastGEAR analysis for the *Bcss + Bt* cluster, the lineages were defined manually to search for potential recombination between *Bcss* and various *Bt* populations.</footnote>

<footnote>** The IDs of the samples that do not belong to any population (i.e., one genome per cluster).</footnote>

<footnote>Note that, as of June 2025, *B. hominis* and *B. weihenstephanensis* have been classified as *B. mycoides* (Parte et al., 2020)</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Maximum likelihood phylogenetic tree of the *B. paranthracis* + *B. pacificus* + *B. tropicus* PopCOGenT cluster 0. The tree was reconstructed using the alignment of core genes and the GTR + F + I + R5 substitution model (Tavaré, 1986). Branch lengths are presented in units of the expected number of nucleotide substitutions per site. Bootstrap values are shown at internal nodes. Each sample label is colored according to the corresponding PopCOGenT subcluster. **A:** PubMLST, predominant sequence types. Jaccard, Groups of genomes with similar accessory gene content. ClonalFrameML, logarithmic sum of recombining regions in the alignment of core genes; Subcluster, PopCOGenT subcluster number. **B:** Admixture blocks (ancestral recombination), patterns of shared ancestry within homologous regions across the core genome alignment. Blocks of the same color represent homologous regions shared by two lineages. The color of each block is assigned based on the lineage with the higher number of samples. FastGEAR detected four lineages: red (*B. pacificus* subclusters 0.1 + 0.2), dark blue (*B. tropicus* subcluster 0.8), green (*B. paranthracis* subclusters 0.6 + 0.7), and blue (*B. paranthracis* subclusters 0.0 + 0.3 + 0.4 + 0.5). White blocks represent recent recombination, indicating putative recombination events occurring in some, but not all, members of a population. **C:** Admixture blocks (recent recombination). The dark blue arrow indicates the sample 11CEB28. The green arrow indicates the sample 18SBCL174. Black blocks indicate putative external recombination. The displayed core genome alignment covers 1.35 Mbp out of the total 2.70 Mbp. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

5.3) and 5.4, spanning from 0.509 to 0.750 Mbp and 1.403 to 1.685 Mbp. Additionally, other admixture regions (0.081–0.114 Mbp, 0.658 to 0.828 Mbp, and 1.033 to 1.073 Mbp) were shared between the (5.1 + 5.2) and (5.0 + 5.3) fastGEAR lineages. Recombination analysis performed by ClonalFrameML, revealed the highest sum of recombining regions in the sample 11CEB61 (0.77 Mbp), followed by 17SBCL1228 (0.61 Mbp), 08CEB128 (0.58 Mbp), 08CEB132BAC (0.55 Mbp), and 09CEB40 (0.52 Mbp) (Fig. 4A). The R/θ, δv, and r/m parameters were estimated at 0.27, 8.9, and 2.4, respectively.

The analysis of accessory gene content identified significant groups that did not follow the population structure (Fig. 4A, Supplementary Fig. S3). Group 12 included samples from subclusters 5.0 and 5.3. Group 13 consisted of three genomes from subcluster 5.4, along with the samples 14SBCL108 and 17SBCL1228 that belonged to subclusters 5.1 and 5.3, respectively.

### 3.5. Toxin genes

We focused primarily on the following *Bcsl* genes: (i) the *ces* operon, involved in cereulide production; (ii) insecticidal *Bt* toxin genes (*cry*, *cyt*, and *vip*); and (iii) the enterotoxin-encoding genes *nheABC*, *hblCDA*, *hlyII* and *cytK*. The *ces* operon was detected in 21 out of 294 samples, with a predominance in *B. paranthracis* subcluster 0.3 (Fig. 1, Supplementary Table S1). In this subcluster, 18 out of 36 genomes contained the operon. The *ces* operon was also identified in *B. paranthracis* subcluster 0.4 (2 out of 10 genomes) and in one *B. mycoides* genome.

The genes encoding insecticidal toxins (Cry, Cyt, and Vip) were detected in all samples within *Bt* subcluster 1.1, except for the four “non-*Bt*” genomes. These genes were also found in some samples of *Bt* cluster 9 (2/8 samples) and *Bt* cluster 10 (1/2). Additionally, they were detected in *Bcss* subcluster 1.0 (3/28), *B. wiedmannii* + “*B. pretiosus*”
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Maximum likelihood phylogenetic tree of the *Bcss* and *Bt* PopCOGenT cluster 1. The tree was reconstructed using the alignment of core genes and the GTR + F + I + R3 substitution model (Tavaré, 1986). Branch lengths are presented in units of the expected number of nucleotide substitutions per site. Bootstrap values are shown at internal nodes. Each sample label is colored according to the corresponding PopCOGenT subcluster. A: PubMLST, predominant sequence types. Jaccard, groups of genomes with similar accessory gene content. ClonaFrameML, logarithmic sum of recombining regions in the alignment of core genes. Subcluster, PopCOGenT subcluster number. B: Admixture blocks (ancestral recombination), patterns of shared ancestry within homologous regions across the core genome alignment. Blocks of the same color represent homologous regions shared by two lineages. The color of each block is assigned based on the lineage with the higher number of samples. FastGEAR detected five lineages: purple (*Bt*), dark blue ("non-*Bt*"), green (*Btk*), red (*Bta*), and blue (*Bcss*). Black blocks indicate putative external recombination. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Maximum likelihood phylogenetic tree of *B. wiedmannii* and "*B. pretiosus*" PopCOGenT cluster 5. The tree was reconstructed using the alignment of core genes and the GTR + F + I + R7 substitution model (Tavaré, 1986). Branch lengths are presented in units of the expected number of nucleotide substitutions per site. Bootstrap values are shown at internal nodes. Each sample label is colored according to the corresponding PopCOGenT subcluster. A: Jaccard, groups of genomes with similar accessory gene content. ClonaFrameML, logarithmic sum of recombining regions in the alignment of core genes. Subcluster, PopCOGenT subcluster number. B: Admixture blocks (ancestral recombination), patterns of shared ancestry within homologous regions across the core genome alignment. Blocks of the same color represent homologous regions shared by two lineages. The color of each block is assigned based on the lineage with the higher number of samples. FastGEAR detected three lineages: green (PopCOGenT subcluster 5.4), red (subclusters 5.0 + 5.3), and blue (subclusters 5.1 + 5.2). White blocks represent recent recombination, indicating putative recombination events occurring in some, but not all, members of a population. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
cluster 5 (3/24), *B. toyonensis* cluster 8 (2/9), and a single *B. albus* genome.

The *cytK-1* gene variant, specific to *B. cytotoxicus* (Fagerlund et al., 2004), was detected in all members of the corresponding cluster 2. The *cytK-2* variant was detected in 6 out of 14 PopCOGenT clusters, mainly occurred in *B. pacificus* subcluster 0.1 (4/4), *B. paranthracis* subcluster 0.4 (9/10), *B. paranthracis* subcluster 0.7 (2/2), *Bcss* + *Bt* cluster 1 (112/113), *B. anthracis* cluster 4 (8/11), *Bt* cluster 9 (4/8), and *Bt* cluster 10 (2/2). The *hlyII* gene was mainly found in *B. wiedmannii* + “B. pretiosus” cluster 5 (21/24), *B. anthracis* cluster 4 (9/11) and *Bt* cluster 9 (7/8).

The *hbl* operon predominantly occurred in *Bcss* + *Bt* cluster 1 (112/113), as well as in almost all members of *B. mycoides* cluster 6 (10/11), *B. toyonensis* cluster 8 (9/9), *Bt* cluster 9 (8/8), and *B. anthracis* cluster 11 (2/2). In *B. wiedmannii* + “*B. pretiosus*” cluster 5, the operon was present in all subclusters, except subcluster 5.4. Lastly, all samples, except 15SBCL1065 and 18SBCL773, carried the *nhe* operon.

### 3.6. Antimicrobial and heavy metal resistance genes

Thereafter, we searched for the genes that might potentially be involved in antibiotic and heavy metal resistance in the Bcsl group. Using AMRFinderPlus, we detected beta-lactamase encoding genes (*bla*1, *bla*2, or *bla*III) in all clusters, except for *B. cytotoxicus* cluster 2 and the unclustered genome 09CEB56 (Fig. 1).

The macrolide resistance genes, *mphL* and *mphM*, encoding macrolide 2'-phosphotransferase, were found in 17.0 % (50/294) and 1.7 % (5/294) of examined samples, respectively. The *mphL* gene was mainly detected in *B. paranthracis* subcluster 0.4 (8/10), *B. anthracis* cluster 4 (11/11), *B. luti* cluster 7 (3/3), *B. toyonensis* cluster 8 (9/9), and *Bt* cluster 9 (5/8). The *tet(45)* gene encoding a tetracycline efflux pump was found in 4.1 % (12/294) of samples, while *tet(L)* was identified in one *B. anthracis* sample (10CEB21). The *vanR* and *vanS* genes, constituting the regulatory component of the *van* operon, were detected in 9.2 % (27/294) of samples. These genes were annotated as the ‘vancomycin resistance response regulator transcription factor VanR’ and the ‘vancomycin resistance histidine kinase VanS,’ respectively. Additionally, we identified the *vanY* and *vanZ* genes, which represented the accessory components of the operon. The *vanY* gene was annotated as ‘VanY domain-containing protein,’ ‘LD-carboxypeptidase,’ or ‘D-alanyl-D-alanine carboxypeptidase’. Notably, the core resistance elements, encoded by the *vanHAX* gene cluster, were absent in all samples. Furthermore, we observed several structural variations (Fig. 5): *vanRSY*, *vanRSYZ*, *vanRS-lipoprotein-vanYZ*, and *vanRSYZ*. In the genome 22SBCL0142_2 (*B. paranthracis* subcluster 0.4), we identified three variants of *vanRSY*. In the first variant, the genes were annotated as vancomycin resistance genes. In the second variant, *vanR* was annotated

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Structural variations in the gene cluster forming regulatory and accessory parts of the putative vancomycin resistance operon detected in *Bcsl*. A: *vanRSY* gene cluster. B: *vanRSYZ* gene cluster. C: *vanRSYZ* gene cluster including a lipoprotein-encoding gene. D: Gene cluster with duplicated *vanY*.</figcaption></figure>

as ‘VanR-ABDEGLN family response regulator transcription factor’, and *vanS* as ‘vancomycin resistance histidine kinase VanS’, with a downstream carboxypeptidase was by the *dacC* gene. In the third variant, *vanR* was annotated as ‘DNA-binding response regulator OmpR family’, *vanS* as ‘Sensor histidine kinase WalK’, and *vanY* as ‘LD-carboxypeptidase LdcB LAS superfamily’.

The mercury (*merB*<sub>3</sub>*R*<sub>1</sub>*EPAR*<sub>2</sub>*B*<sub>2</sub>*B*<sub>1</sub>) (Naguib et al., 2018) and arsenic (*arsRBCDA*) (Chien et al., 2019) resistance operons were observed in 3.1 % (9/294) and 17.0 % (50/294) of genomes, respectively (Fig. 1).

Finally, we reviewed the type of dishes associated with strains carrying antibiotic resistance genes (Supplementary Table S1). We noticed that the isolates carrying the *tet* genes were much more frequently associated with dishes marinated and/or containing spices (62 % versus 7 %) than the *tet*-negative isolates. Similarly, although with a smaller gap, we observed that the samples carrying *mphL* were more frequently associated with the “Fish, shellfish, and their products” FAO Level 1 category (16 % versus 5 %).

### 3.7. Antibiotic susceptibility phenotype

To confirm the phenotype of loss of antibiotic susceptibility due to the presence of previously highlighted antibiotic resistance genes (ARGs), we carried out disk diffusion assays on 64 Bcsl isolates (Fig. 6, Supplementary Table S4). We first compared the susceptibility to ampicillin (a β-lactam class antibiotic) between *bla*+ (n = 53) and *bla*- (n = 11) isolates. The isolates not carrying *bla* genes, predominantly belonging to *B. cytotoxicus* (10/11), demonstrated a large and significant increase in ampicillin susceptibility compared to *bla*+ isolates. Concerning erythromycin, a macrolide antibiotic, results revealed that the isolates carrying *mphL* (n = 17) exhibited a significantly reduced inhibition diameter, compared to the *mphL*-isolates (n = 47). According to the breakpoint defined by EUCAST (R < 24 mm), 88 % (15/17) of *mphL*+ isolates were resistant to macrolides. Similarly, Bcsl isolates carrying the *tet(45)* gene (n = 12) exhibited a significantly reduced inhibition diameter when being exposed to tetracyclin (a cyclin class antibiotic), compared to the *tet(45)*- isolates (n = 51). Moreover, a single isolate carrying the *tet(L)* gene showed a total lack of sensitivity to this antibiotic. Finally, the isolates carrying the *vanR* and *vanS* genes (n = 20) showed a significant, albeit small, reduction in inhibition diameter than the isolates lacking these genes (n = 44), after exposure to vancomycin, a glycopeptide-belonging antibiotic. However, according to EUCAST criteria (R < 10 mm), all the isolates were susceptible to vancomycin.

### 3.8. Proportions of Bcsl populations

To search for associations between the most frequently found *Bcsl* populations and food categories, we applied a bootstrap-based approach to reflect the true proportions of the FBO-associated *Bcsl* populations in France. The most frequently detected populations were *Bcss* subcluster 1.0 (17.0 ± 0.3 %), *B. paranthracis* subcluster 0.3 (16.1 ± 1.0 %), and *Btk* subcluster 1.1 (7.6 ± 1.0 %) (Fig. 7). In FBOs where a single isolate was detected, the most frequent population was *B. paranthracis* subcluster 0.3, accounting for 6.3 ± 1.4 % of all FBO cases, followed by *Btk* (3.0 ± 1.3 %) and *Bcss* (2.8 ± 0.5 %).

### 3.9. Distribution of Bcsl populations across food categories

Thereafter, we estimated the distribution of the entire *Bcsl* group and the three predominant populations across food categories. Statistically significant associations were found when using FoodEx2 classification (European Food Safety Authority, 2015) (Supplementary Table S1). At Level 1, the *Bcsl* group was most prevalent in the “Composite dishes” category (51 ± 0.1 %), followed by “Grains and grain-based products” (14 ± 0.1 %), and “Vegetables and vegetable products” (14 ± 0.1 %) (Fig. 8A). At Level 2, *Bcsl* was primarily found in “Dishes including ready-to-eat meals excluding soups and salads” (39 ± 0.1 %) and “Soups
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Antibiotic susceptibility testing of *Bcsl* performed using the Kirby-Bauer disk diffusion method (EUCAST method v13.0 and breakpoint tables for interpretation v15.0). For each antibiotic, isolates were selected based on the presence or absence of the identified antimicrobial resistance genes of interest. A: Resistance to ampicillin (10 µg, AMP10) examined in the *bla*- and *bla* + isolates; *bla*, beta-lactamase encoding genes. B: Resistance to erythromycin (15 µg, E15) resistance examined in the *mphL*- and *mphL* + isolates; *mphL*, macrolide 2'-phosphotransferase encoding gene. C: Resistance to tetracycline (30 µg, TE30) examined for the *tet*-, *tet(45)*+ and *tet(L)* + isolates; *tet*, tetracycline efflux pump encoding gene. D: Resistance to vancomycin (5 µg, VA5) examined in the *van*- and *van* + isolates; *van*-putative regulatory part (*vanRS*) of the vancomycin resistance operon. For B and D, dotted lines indicate the respective values of EUCAST zone diameter breakpoints (Breakpoint Table v.14.0). ****, *p*-value<0.0001 (Mann-Whitney test). The figure was designed using GraphPad v.10.4.1.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Proportions of *Bcsl* populations estimated from the collection of sequenced genomes sampled during foodborne outbreaks (FBO) in France between 2004 and 2023. **Presumptive Bcsl**, the sample is counted even if another pathogen (another species or mix of *Bcsl* species) is identified in the FBO. **One Bcsl**, only a single *Bcsl* isolate detected in the FBO. *Bcss*, *B. cereus sensu stricto* subcluster 1.0. *Bpar* 0.3, *B. paranthracis* subcluster 0.3. *Btk*, *Bt* subsp. *kurstaki* subcluster 1.1. 4, *B. anthracis* cluster 4. 8, *B. toyonensis* cluster 8. 0.4, *B. paranthracis* subcluster 0.4. 5.1, *B. wiedmannii* + "*B. pretiosus*" subclusters 5.1 + 5.2. 6, *B. mycoides* cluster 6. Unclust., unclustered genome. 9, *Bt* cluster 9. 0.2, *B. pacificus* subcluster 0.2. 5.0, *B. wiedmannii* + "*B. pretiosus*" subclusters 5.0 + 5.3 + 5.4. *Bta*, *Bt* subsp. *aizawai* subcluster 1.1. 0.0, *B. paranthracis* subclusters 0.0 + 0.5. 1.1 non-*Bt*, particular *Bt* strains from subcluster 1.1. 2, *B. cytotoxicus* cluster 2. 0.1, *B. pacificus* subcluster 0.1. 0.6, *B. paranthracis* subcluster 0.6. 3, *B. mobilis* cluster 3. 0.8, *B. tropicus* subcluster 0.8. 7, *B. luti* cluster 7. 11, *B. anthracis* cluster 11. 0.7, *B. paranthracis* subcluster 0.7. 12, *B. paramobilis* cluster 12; 13, *B. wiedmannii* cluster 13. 10, *Bt* cluster 10. 1.1 *Bt*, *Bt* ST13 from subcluster 1.1.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Distribution of the entire *Bcsl* group and the predominant populations across different food categories determined using FoodEx2 exposure classification Level 1 and 2 (European Food Safety Authority, 2015). **A: FoodEx Level 1.** Composite_dishes, "Composite dishes" category. Fish_seafood, "Fish, seafood, amphibians, reptiles and invertebrates". Fruit, "Fruit and fruit products". Grains, "Grains and grain-based products". Legumes_nuts_spices, "Legumes, nuts, oilseeds and spices". Meat, "Meat and meat products". Milk, "Milk and dairy products". Seasoning_spices, "seasoning sauces and condiments". Starchy_rootsl, "Starchy roots or tubers and products thereof sugar plants". Sugar, "Sugar and similar confectionery and water-based sweet desserts". Vegetables, "Vegetables and their products". **B: FoodEx Level 2.** Cereals, "Cereals and cereal primary derivatives". Cheese, "Cheese". Crustaceans, "Crustaceans". Dishes_excl_soups_plads, "Dishes including ready-to-eat meals excluding soups and salads". Fine_bakery, "Fine bakery wares". Fish, "Fish meat". Flower_brassica, "Flowering brassica". Fungi, "Fungi, mosses andlichens". Leafy_vegetables, "Leafy vegetables". Milk, "Milk and dairy powders and concentrates". Mollusks, "Mollusks". Processed_fruit, "Processed fruit products". Processed_vegetables, "Processed or preserved vegetables and similar". Root_excl_starch, "Root and tuber vegetables excluding starchy and sugar". Sausages, "Sausages". Soups_plads, "Soups and salads". Spices, "Spices". Starchy_rootsl, "Starchy roots and tubers". Sugar_plants, "Sugar plants". *Bcsl*, the *B. cereus sensu lato* entire group. *Bcss*, *B. cereus sensu stricto* subcluster 1.0. *Bpar*, *B. paranthracis* subcluster 0.3. *Btk*, *Bt* subsp. *kurstaki* subcluster 1.1.</figcaption></figure>

and salads" (15 ± 0.1 %), both subcategories of "Composite dishes", followed by "Cereals and cereal primary derivatives" (13 ± 0.1 %), a subcategory of "Grains and grain-based products" (Fig. 8B). At Level 1, analysis of the three predominant populations revealed a significant positive association between *Bcss* and "Composite dishes" (*p* = 0.029), where it was present in 73 ± 3.9 % of cases (Fig. 8A). *B. paranthracis* subcluster 0.3 was detected in "Grains and grain-based products" and "Starchy roots and tubers" in 26 ± 5.8 % and 10 ± 3.7 % of cases, respectively, though these associations were not statistically significant when comparing with the overall *Bcsl* background. Similarly, *Btk* was found in "Fruit and fruit products" in 13 ± 6.4 % of cases, but without statistical significance. At Level 2, *B. paranthracis* subcluster 0.3 showed a significant positive association with "Cereals and cereal primary derivatives" (*p* = 0.033), accounting for 24 ± 5.5 % of cases. Additionally, *Btk* was significantly associated with "Soups and salads" (34 ± 12.2 %) (*p* < 0.0001), which primarily included raw vegetable-based salads (Fig. 8B, Supplementary Fig. S1).

4. Discussion

In this study, we resolved the population structure of *Bcsl* collected
during FBOs in France between 2004 and 2023.

First, we observed putative signatures of horizontal transfer for the *cry/cyt/vip* genes and the *ces* operon. Indeed, horizontal transfer is well-documented for *Bcsl*. For instance, although plasmid-associated genes encoding insecticidal toxins (Cry, Cyt, and Vip) are characteristic of *Bt* (*Schnepf et al., 1998*), we also identified them in the other species, comprising three cases of *Bcss*. Consistently, other *Bcsl* species possessing insecticidal toxin genes were described in various studies (Lazarte et al., 2018; Sauka et al., 2022). Sauka et al. (2022) demonstrated that the *B. toyonensis* strain Bto-UNVM_94, which carries a *cry* gene, can produce parasporal crystals with insecticidal activity. Similarly, the plasmid-located *ces* operon (*Ehling-Schulz et al., 2006*), encoding genes involved in the synthesis of the emetic toxin, is mostly found in a part of *B. paranthracis* ST26 samples, as also revealed in previous studies (Carroll et al., 2022b; Frentzel et al., 2024). We revealed one case of *B. mycoides* carrying the *ces* operon. Similarly, several *B. mycoides* strains
isolated from soil were shown to carry the *cesHPTABCD* operon as part of the composite transposon *Tnces* (Mei et al., 2014). *Frentzel et al.* (2024) described an emetic *B. mobilis* strain isolated from meat salad, able to produce cereulide. In the study of *Pheepakpraw et al.* (2023), the presence of the *cesA* and *cesB* gene fragments was confirmed for several FBO- and clinic-associated strains, including *Bt* ST8, *Bt* ST15, and *Bcsl* ST2804 and ST2805. For the other toxin-encoding genes, we demonstrated that *nhe* was present in almost all samples, while *hbl*, *hlyII* and *cytK-2* were found in various lineages. This supports the findings of Böhm et al. (2015), who showed that the *nhe* operon was present in all *Bcsl* and was vertically inherited with extremely rare cases of gene loss. In contrast, *hbl* and *cytK-2* were shown to be subject to horizontal transfer (Böhm et al., 2015). Finally, we observed the predominant occurrence of the *hbl* operon in *panC* phylogenetic groups II, IV, V, and VI. This aligns with the findings of Carroll et al. (2017), who demonstrated a statistically significant association between the presence of the operon and these clades.

Furthermore, we sporadically identified the antimicrobial and heavy metal resistance genes, such as *tet(45)*, *tet(L)*, *mphL*, *ars*, and *mer*, within distinct *Bcsl* populations, suggesting the adaptation of *Bcsl* to specific environmental conditions. We noticed that the *Bcsl* isolates carrying the *tet* genes were much more frequently detected in dishes containing spices. Given that most spices in France are imported, we speculate that this may reflect a potential spread of ARG-carrying *Bcsl* across geographical borders, particularly from regions with higher antibiotic pressure (Okeke and Edelman, 2001). Similarly, the *mph* genes were frequently associated with seafood products, for which the spread of ARGs is a well-documented concern, both in marine environments (Ferraro et al., 2024) and fish farming areas (Helsens et al., 2020).

By contrast, we found that all *Bcsl*, except *B. cytotoxicus* cluster 2, possessed genes encoding various beta-lactamases. Indeed, numerous studies demonstrated 98–100 % resistance of *Bcsl* isolates to penicillin, cefotaxime, ampicillin and oxacillin, concluding natural resistance of *Bcsl* to beta-lactams (Glasset et al., 2018; Bianco et al., 2021; Jung et al., 2022; Rajalingam et al., 2022; Mohammadi et al., 2023). To date, this is the first time that a specific *Bcsl* population (*B. cytotoxicus*) is exhaustively associated with the absence of beta-lactamases and shows a drastically increased sensitivity phenotype to a beta-lactam antibiotic, compared to the rest of the group.

Thereafter, when searching for the *van* operon, we detected several variants of the *vanRSYZ* gene cluster with no resistance-associated genes (*vanH* encoding dehydrogenase, *vanA/B/D/F/M* encoding D-ala-D-lactate ligase, and *vanX* encoding D,D-dipeptidase) (Stogios and Savchenko, 2020). Indeed, the phylogenetic trees reconstructed by Kardos et al. (2024) for the vancomycin resistance genes in Bacilli and Clostridia sequences available in the NCBI database showed that the acquisition of *vanRS*, with or without *vanYZ*, could occur in several clades. These acquisitions may either be associated with resistance-related genes or exist independently, suggesting that *vanRS*-like regulatory systems could potentially regulate other genes unrelated to the *van* operon (Kardos et al., 2024). This is also supported by numerous studies showing 95–100 % susceptibility of *Bcsl* strains to vancomycin (Glasset et al., 2018; Bianco et al., 2021; Jung et al., 2022; Rajalingam et al., 2022; Mohammadi et al., 2023).

To sum up, future research should complement gene presence-absence analysis with gene expression analysis, genetic and structural variation analysis, and genome-wide association studies to link specific genotypes with phenotypic traits, such as toxin production levels and antibiotic resistance. This is crucial for the *Bcsl* group in the context of public health, as non-pathogenic strains may acquire virulence factors and develop antibiotic resistance via horizontal gene transfer and mutation.

The clustering of *Bcsl* strains enabled us to detect three composite PopCOGenT clusters, each representing a group of closely related species: *B. paranthracis* + *B. pacificus* + *B. tropicus* cluster 0, *Bcsl* + *Bt* cluster 1, and *B. wiedmannii* + "*B. pretiosus*" cluster 5. In these clusters,

FastGEAR identified numerous admixture blocks within lineages, suggesting either shared ancestry or recombination. Metrics from Clonal-FrameML show that recombination has a greater overall impact than mutation in the examined *Bcsl* clusters, representing an intermediate level between clonal and true recombining bacterial species (Vos and Didelot, 2008; Patino-Navarrete and Sanchis, 2017). This suggests that recombination may serve as an important mechanism for *Bcsl* to share fitness-enhancing alleles.

Regarding cluster 5, the *B. wiedmannii* type strain was originally isolated from raw milk and characterized as a psychrotolerant isolate producing hemolysin BL and non-hemolytic enterotoxin, as well as showing cytotoxicity to HeLa cells (Miller et al., 2016). "*B. pretiosus*", effectively published as species phylogenetically close to *B. wiedmannii*, was initially isolated from the rhizosphere of alfalfa (Robas Mora et al., 2023). We demonstrated that significant Groups based on the Jaccard index did not entirely align with the population structure. This suggests that the *B. wiedmannii* + "*B. pretiosus*" cluster occupies a broad ecological niche, with strains adapting to local environmental conditions through a unique set of accessory genes.

In *Bcsl* + *Bt* cluster 1, *Bt* subcluster 1.1 comprised several *Bt* lineages: *Bta*, *Btk*, *Bt* ST13, and the group of samples lacking *Bt*-specific genetic markers, as well as *cry*, *cyt*, or *vip* genes. Nevertheless, these samples are concluded to be particular *Bt* strains, as demonstrated by phylogenetic and admixture analyses. In addition, we want to emphasize that the study by Fichant et al. (2022) aimed to develop the workflow specifically for pesticide and FBO-associated *Bta* and *Btk* strains rather than for the full diversity of *Bt*. The *Bcsl* subcluster 1.0 is represented by a single PopCOGenT/fastGEAR population showing complex population structure. The subcluster includes various STs, with predominant ST24. In *Bcsl* ST24 and *Bta*, specific sets of accessory genes were enriched in the "Mobilome: phages, transposons" functional category. Indeed, the mobilome of *Bcsl* includes various transposable elements (TEs), such as insertion sequences and class II transposons that are often associated with virulence genes (Mei et al., 2014; Fayad et al., 2019).

Finally, we found that *Bt* subcluster 1.1 was distinct from *Bt* clusters 9 and 10, thereby showing *Bt* to be a polyphyletic group. The *Bt* samples from clusters 9 and 10 represent various *Bt* sequence types, including *Bt* serovar *israelensis* (ST16 and ST197) and *Bt* serovar *morrisoni* (ST23) (Supplementary Table S1) (Wang et al., 2018). Notably, only three out of ten strains carry delta-endotoxin-encoding genes. This observation suggests the loss or sporadic acquisition of insecticidal toxin-encoding genes within these *Bt* phylogenetic groups.

Therefore, genetic and phenotypic overlap within the *Bcsl* + *Bt* cluster may complicate the differentiation between *Bt* and *Bcsl* strains carrying *cry* genes, as well as between environmental *Bt* strains and pesticidal strains, using laboratory methods such as PCR of the *panC* gene or microscopic observation of insecticidal crystals. An accurate identification of pesticidal *Bt* remains crucial, given the increasing evidence that *Bt* contaminates food products (EFSA Panel on Biological Hazards, 2016; Frentzel et al., 2024; Carroll et al., 2022b), disseminates in hospital environments (Muigg et al., 2022), and is associated with FBOs (Bonis et al., 2021; Biggel et al., 2022a, 2022b). Nevertheless, the initial steps to detect pesticidal *Bt* strains, primarily derived from *Bta* and *Btk* subspecies, were performed by Frentzel et al. (2020), Bonis et al. (2021), Biggel et al. (2022a), and Fichant et al. (2022). Frentzel et al. (2020) found that some *Bt* strains isolated from vegetables differed from pesticidal strains by ≤ 4 wgSNPs. Biggel et al. (2022a) demonstrated that the genetic distance between several food- and FBO-associated *Bta* and *Btk* strains and pesticidal strains was no more than 2 cgSNPs or 6 wgSNPs, thereby suggesting a pesticidal origin for the food- and FBO-associated strains. Bonis et al. (2021) showed that pesticidal and FBO-associated *Bt* strains formed four phylogenetic clusters (a to d), with the median cgSNP distances ranging from 0 to 2 within each cluster. However, to date, there is no accepted threshold of genetic distance between pesticidal and non-pesticidal *Bt* genomes to classify them as the same strain. Subsequently, Fichant et al. (2022) proposed
the approach to detect accessory genes specific for pesticidal and FBO-associated *Bt* strains at the species level (*Bt* and non-*Bt*), subspecies level (*Bta* and *Btk*), and cluster level (a to d). Further precise detection of *Bt* strains remains a subject for future work, requiring the identification of finer genetic markers such as SNPs or indels.

The largest cluster, *B. paranthracis* + *B. pacificus* + *B. tropicus* cluster 0, comprises three closely related species. Among them, *B. paranthracis* subcluster 0.3 contains the highest number of samples and consists of several strains, with predominant ST26 and ST266. Complex population structure of this cluster suggests a broad ecological niche, allowing the populations to adapt to a wide range of environments. The three species were initially isolated from ocean sediment (Liu et al., 2017). *B. paranthracis* is often associated with foodborne outbreaks (Messelhäusser et al., 2014). *B. pacificus* is found in soil, with some strains showing plant-growth promoting properties (Ma et al., 2024). *B. tropicus* is also found in soil and sludge environments (Thakur et al., 2021; Shen et al., 2022). Recombination likely plays a significant role in shaping the population structure of this cluster, especially since the R/θ and r/m ratios were higher compared to those for clusters 1 and 5. We revealed several strains, such as 11CEB28 (*B. tropicus* subcluster 0.8), 18SBCL174 (*B. paranthracis* subcluster 0.7) and 22SBCL0891.1 (*B. paranthracis* subcluster 0.7), that demonstrated the highest number of recent recombination blocks and the longest recombining regions, suggesting they are products of recent recombination between distinct populations or have recently inherited these signatures.

The analysis of accessory gene content revealed two distinct groups, Group 33 and Group 41, within *B. paranthracis* subcluster 0.3, each possessing specific accessory genes likely suited to specific environmental conditions. Interestingly, the accessory genes of *B. paranthracis* subclusters 0.0 and 0.5, and the clinical samples GCF_021014465, GCA_029696605, and GCF_036851755, are enriched in 'Cell wall/membrane/envelope biogenesis' COG functional category. The genes from this functional category are known to play a role in the synthesis of cell wall components, including virulence factors, and to contribute to host-pathogen interactions (King and Roberts, 2016).

The detailed population structure analysis helped us accurately differentiate *Bcsl* samples into distinct populations and estimate the prevalence of each population across different food categories. The most frequent FBO-associated populations were *Bcss*, *B. paranthracis* subcluster 0.3 including emetic strains, and *Btk*. The predominance of the commonly known foodborne pathogens *Bcss* and *B. paranthracis* was also demonstrated in other studies (Messelhäusser et al., 2014; Frentzel et al., 2024; Jacobs et al., 2025). The predominance of *Btk* can be explained by the widespread use of *Bt* pesticide strains and subsequent contamination of food, particularly vegetable-based meals, which aligns with the findings of Bonis et al. (2021). Furthermore, *B. paranthracis* subcluster 0.3 showed a significant positive association with starchy foods, such as cereals, while *Bcss* was associated with composite dishes. These findings align with the study by Messelhäusser et al. (2014), who reported that cereulide-producing strains isolated during FBOs in Bavaria between 2007 and 2013 predominantly occurred in cereals. Additionally, various reports (Osimani et al., 2018; Santé Publique France. Bulletin, 2022) highlighted that the *Bcsl* group was commonly linked to composite dishes, including pasta- and rice-based dishes. This indicates the elevated risks associated with these food categories and highlights the challenge of tracing the initial source of contamination.

## 5. Conclusion

In our study, we resolved the population structure of the large collection of FBO-associated *Bcsl* strains collected in France between 2004 and 2023. Through the analysis of 294 genomes from 183 FBOs, we identified 14 PopCOGenT clusters corresponding to validly and effectively published species. Three of these clusters exhibited a composite structure. The predominant FBO-associated populations in France are *Bcss* subcluster 1.0 including various STs, *B. paranthracis* subcluster

0.3 including emetic strains, and *Btk* subcluster 1.1. Each population was positively associated with the specific food category. We believe that the reverse ecology approach applied in this study proved effective to accurately differentiate *Bcsl* strains and obtain biologically meaningful clusters and subclusters. We consider this method to be suitable for future epidemiological studies. It allows for statistically meaningful sample sizes, helps clarify the confusion around *Bcsl* taxonomy, and avoids the use of arbitrary thresholds in genomospecies definition. The predominant populations identified in this study will be used in further dose-response modeling and will be associated with epidemiological data, such as the percentage of affected individuals, age, symptoms, infectious dose, and collection point. Furthermore, other key factors, such as the presence-absence of toxin-encoding genes, toxin production levels, and cytotoxicity, can be incorporated as additional variables for specific populations in risk modeling.

### CRediT authorship contribution statement

**Ksenia Mozhaitseva:** Writing – original draft, Visualization, Investigation, Formal analysis, Conceptualization. **Sylvie Pairaud:** Methodology, Investigation. **Olivier Firmesse:** Writing – review & editing, Validation, Supervision, Conceptualization. **Mathilde Bonis:** Writing – review & editing, Validation, Supervision, Investigation, Funding acquisition, Conceptualization.

### Data availability statement

The accessions for the raw reads of the studied *Bcsl* strains are provided in Supplementary Table S1.

### Funding

The study was performed as part of the National Agency for Research project "Virulence Assessment and Innate Immune Response Monitoring of *Bacillus cereus* Spores During Infection – BaDAss" (ANR-22-CE35-0006).

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

We thank all the district veterinary and food analysis laboratories for carrying out *Bcsl* detection and transmitting isolates together with epidemiologic data to the Food Safety Laboratory. We thank Virginie Chesnais and Déborah Merda (Shared Data Analysis Support Service, ANSES) for developing of the bacflow assembly pipeline and their assistance in the management of software tools. We also thank Sophie Granier for her expertise in antibiotic resistance testing and Laurent Guillier for the fruitful discussions on the statistical and phylogenetic aspects of the project.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.fm.2025.104882.

### References

Arevalo, P., VanInsberghe, D., Elsherbini, J., Gore, J., Polz, M.F., 2019. A reverse ecology approach based on a biological definition of microbial populations. Cell 178 (4), 820–834. https://doi.org/10.1016/j.cell.2019.06.033.
Barakat, S., Kim, H., Dankar, R., Hewlett, C., 2024. *Bacillus thuringiensis* bacteremia in a 30-Year-Old intravenous drug user: a report of a rare case. Cureus 16 (10), e71704. https://doi.org/10.7759/cureus.71704.

Battisti, L., Green, B.D., Thorne, C.B., 1985. Mating system for transfer of plasmids among *Bacillus anthracis*, *Bacillus cereus*, and *Bacillus thuringiensis*. J. Bacteriol. 162 (2), 543–550. https://doi.org/10.1128/jb.162.2.543-550.1985.

Bazinet, A.L., 2017. Pan-genome and phylogeny of *Bacillus cereus* sensu lato. BMC Evol. Biol. 17, 1–16. https://doi.org/10.1186/s12862-017-1020-1.

Bianco, A., Capozzi, L., Monno, M.R., Del Sambro, L., Manzulli, V., Pesole, G., Parisi, A., 2021. Characterization of *Bacillus cereus* group isolates from human bacteremia by whole-genome sequencing. Front. Microbiol. 11. https://doi.org/10.3389/fmicb.2020.599524. Article 599524.

Biggel, M., Etter, D., Corti, S., Brodmann, P., Stephan, R., Ehling-Schulz, M., Johler, S., 2022a. Whole genome sequencing reveals biopesticidal origin of *Bacillus thuringiensis* in foods. Front. Microbiol. 12. https://doi.org/10.3389/fmicb.2021.775669. Article 775669.

Biggel, M., Jessberger, N., Kovac, J., Johler, S., 2022b. Recent paradigm shifts in the perception of the role of *Bacillus thuringiensis* in foodborne disease. Food Microbiol. 105. https://doi.org/10.1016/j.fm.2022.104025. Article 104025.

Böhm, M.E., Huptas, C., Krey, V.M., Scherer, S., 2015. Massive horizontal gene transfer, strictly vertical inheritance and ancient duplications differentially shape the evolution of *Bacillus cereus* enterotoxin operons *hbl*, *cytK* and *nhe*. BMC Evol. Biol. 15, 1–17. https://doi.org/10.1186/s12862-015-0529-4.

Bonis, M., Felten, A., Pairaud, S., Dijoux, A., Maladen, V., Mallet, L., Herbin, S., 2021. Comparative phenotypic, genotypic and genomic analyses of *Bacillus thuringiensis* associated with foodborne outbreaks in France. PLoS One 16 (2). https://doi.org/10.1371/journal.pone.0246885. Article e0246885.

Brynlildsrud, O., Bohlin, J., Scheffer, L., Eldholm, V., 2016. Rapid scoring of genes in microbial pan-genome-wide association studies with Scoary. Genome Biol. 17, 1–9. https://doi.org/10.1186/s13059-016-1108-8.

Carlson, C.R., Johansen, T., Kolsto, A.B., 1996. The chromosome map of *Bacillus thuringiensis* subsp. *canadensis* HD24, is highly similar to that of the *Bacillus cereus* type strain ATCC 14579. FEMS Microbiol. Lett. 141 (2–3), 163–167. https://onlineli brary.wiley.com/doi/pdf/10.1111/j.1574-6968.1996.tb08379.x.

Carroll, L.M., Kovac, J., Miller, R.A., Wiedmann, M., 2017. Rapid, high-throughput identification of anthrax-causing and emetic *Bacillus cereus* group genome assemblies via BTyper, a computational tool for virulence-based classification of *Bacillus cereus* group isolates by using nucleotide sequencing data. Appl. Environ. Microbiol. 83 (17). https://doi.org/10.1128/AEM.010906-17. Article e010906-17.

Carroll, L.M., Cheng, R.A., Kovac, J., 2020a. No assembly required: using BTyper3 to assess the congruency of a proposed taxonomic framework for the *Bacillus cereus* group with historical typing methods. Front. Microbiol. 11. https://doi.org/10.3389/fmicb.2020.580691. Article 580691.

Carroll, L.M., Wiedmann, M., Kovac, J., 2020b. Proposal of a taxonomic nomenclature for the *Bacillus cereus* group which reconciles genomic definitions of bacterial species with clinical and industrial phenotypes. mBio 11 (1). https://doi.org/10.1128/ mbio.00034-20. Article 10-1128.

Carroll, L.M., Cheng, R.A., Wiedmann, M., Kovac, J., 2022a. Keeping up with the *Bacillus cereus* group: taxonomy through the genomics era and beyond. Crit. Rev. Food Sci. Nutr. 62 (28), 7677–7702. https://doi.org/10.1080/10408398.2021.1916735.

Carroll, L.M., Pierneef, R., Mathole, A., Atanda, A., Matle, I., 2022b. Genomic sequencing of *Bacillus cereus* Sensu Lato strains isolated from meat and poultry products in South Africa enables inter- and intranational surveillance and source tracking. Microbiol. Spectra. 10 (3). https://doi.org/10.1128/spectra.000700-22. Article e00700-22.

Chien, M.F., Ho, Y.N., Yang, H.E., Narita, M., Miyauchi, K., Endo, G., Huang, C.C., 2019. Identification of a novel arsenic resistance transposon nested in a mercury resistance transposon of *Bacillus* sp. MB24. Microorganisms 7 (11). https://doi.org/10.3390/ microorganisms7110566. Article 566.

Didelet, X., Wilson, D.J., 2015. ClonalFrameML: efficient inference of recombination in whole bacterial genomes. PLoS Comput. Biol. 11 (2), e1004041. https://doi.org/10.1371/journal.pcbi.1004041.

Edgar, R.C., 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32 (5), 1792–1797. https://doi.org/10.1093/nar/ gkh340.

EFSA Panel on Biological Hazards (BIOHAZ), 2016. Risks for public health related to the presence of *Bacillus cereus* and other *Bacillus* spp. including *Bacillus thuringiensis* in foodstuffs. EFSA J. 14 (7). https://doi.org/10.2903/j.efsa.2016.4524. Article e04524.

Ehling-Schulz, M., Fricker, M., Grallert, H., Rieck, P., Wagner, M., Scherer, S., 2006. Cereulide synthetase gene cluster from emetic *Bacillus cereus*: structure and location on a mega virulence plasmid related to *Bacillus anthracis* toxin plasmid pX01. BMC Microbiol. 6, 1–11. https://doi.org/10.1186/1471-2180-6-20.

Ehling-Schulz, M., Lereclus, D., Koehler, T.M., 2019. The *Bacillus cereus* group: *bacillus* species with pathogenic potential. Microbiol. Spectra. 7 (3). https://doi.org/10.1128/ microbiolspec.gpp3-0032-2018. Article 10-1128.

Ehling-Schulz, M., Knutsson, R., Scherer, S., 2010. *Bacillus cereus*. Genomes of Foodborne and Waterborne Pathogens, pp. 147–164. https://doi.org/10.1128/ 9781555816902.ch11.

Emms, D.M., Kelly, S., 2015. OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthologous inference accuracy. Genome Biol. 16, 1–14. https://doi.org/10.1186/13059-015-0721-2.

European Food Safety Authority (EFSA), 2015. The food classification and description system FoodEx 2 (revision 2). EFSA J. 12 (5), 1–90. https://efsa.onlinelibrary.wiley. com/doi/epdf/10.2903/sp.efsa.2015.EN-804.

European Food Safety Authority (EFSA), European Centre for Disease Prevention and Control (ECDC), 2024. The European Union health 2023 zoonoses report. EFSA J. 22 (12). https://doi.org/10.2903/j.efsa.2024.9106. Article e9106.

Fagerlund, A., Ween, O., Lund, T., Hardy, S.P., Granum, P.E., 2004. Genetic and functional analysis of the *cytK* family of genes in *Bacillus cereus*. Microbiology 150 (8), 2689–2697. https://doi.org/10.1099/mic.0.26975-0.

Fayad, N., Kallassy Awad, M., Mahillon, J., 2019. Diversity of *Bacillus cereus* sensu lato mobilome. BMC Genom. 20, 1–11. https://doi.org/10.1186/s12864-019-5764-4.

Feldgarden, M., Brover, V., Gonzalez-Escalona, N., Frye, J.G., Haendiges, J., Haft, D.H., et al., 2021. AMRFinderPlus and the Reference gene catalog facilitate examination of the genomic links among antimicrobial resistance, stress response, and virulence. Sci. Rep. 11. https://doi.org/10.1038/s41598-021-91456-0. Article 12728.

Ferraro, G.B., Brandtner, D., Franco, A., Iaconelli, M., Mancini, P., Veneri, C., et al., 2024. Global quantification and distribution of antibiotic resistance genes in oceans and seas: anthropogenic impacts and regional variability. Sci. Total Environ. https://doi. org/10.1016/j.scitotenv.2024.176765. Article 176765.

Fichant, A., Felten, A., Gallet, A., Firmesse, O., Bonis, M., 2022. Identification of genetic markers for the detection of *Bacillus thuringiensis* strains of interest for food safety. Foods 11 (23). https://doi.org/10.3390/foods11233924. Article 3924.

Fichant, A., Lanceleur, R., Hachfi, S., Brun-Barale, A., Blier, A.L., Firmesse, O., et al., 2024. New approach methods to assess the enteropathogenic potential of strains of the *Bacillus cereus* group, including *Bacillus thuringiensis*. Foods 13 (8). https://doi. org/10.3390/foods13081140. Article 1140.

Flügge, C., 1886. Die Mikroorganismen. Leipzig: Vogel.

Frentzel, H., Juraschek, K., Pauly, N., Kelner-Burgos, Y., Wichmann-Schauer, H., 2020. Indications of biopesticidal *Bacillus thuringiensis* strains in bell pepper and tomato. Int. J. Food Microbiol. 321. https://doi.org/10.1016/j.ijfoodmicro.2020.108542. Article 108542.

Frentzel, H., Kraemer, M., Kelner-Burgos, Y., Velze, L., Bodi, D., 2024. Cereulide production capacities and genetic properties of 31 emetic *Bacillus cereus* group strains. Int. J. Food Microbiol. 417. https://doi.org/10.1016/j. ijfoodmicro.2024.110694. Article 110694.

Glasset, B., Herbin, S., Granier, S.A., Cavalié, L., Lafeuille, E., Guérin, C., et al., 2018. *Bacillus cereus*, a serious cause of nosocomial infections: epidemiologic and genetic survey. PLoS One 13 (5). https://doi.org/10.1371/journal.pone.0194346. Article e0194346.

Guinebretière, M.H., Thompson, F.L., Sorokin, A., Normand, P., Davyndt, P., Ehling-Schulz, M., et al., 2008. Ecological diversification in the *Bacillus cereus* group. Environ. Microbiol. 10 (4), 851–865. https://doi.org/10.1111/j.1462- 2920.2007.01495.x.

Guinebretière, M.H., Auger, S., Galleron, N., Contzen, M., De Sarrau, B., De Buyser, M.L., et al., 2013. *Bacillus cytotoxicus* sp. Nov. is a novel thermotolerant species of the *Bacillus cereus* group occasionally associated with food poisoning. Int. J. Syst. Evol. Microbiol. 63 (Pt 1), 31–40. https://doi.org/10.1099/ijss.0.0306270. Article e0194346.

Helsens, N., Calvez, S., Prevost, H., Bouju-Albert, A., Maillet, A., Rossero, A., et al., 2020. Antibiotic resistance genes and bacterial communities of farmed rainbow trout fillets (*Oncorhynchus mykiss*). Front. Microbiol. 11. https://doi.org/10.3389/ fmicb.2020.590902. Article e59092.

Hoang, D.T., Chernomor, O., Von Haeseler, A., Minh, B.Q., Vinh, L.S., 2018. UFBoot2: improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35 (2), 518–522. https://doi.org/10.1093/molbev/msx281.

Jacobs, B., Bogaerts, B., Verhaegen, M., Vanneste, K., De Keersmaecker, S.C., Roosens, N. H., Van Hoorde, K., 2025. Whole-genome sequencing of soil- and foodborne *Bacillus cereus* sensu lato indicates no clear association between their virulence repertoire, genomic diversity and food matrix. Int. J. Food Microbiol. 439. https://doi.org/ 10.1016/j.ijfoodmicro.2025.111266. Article 111266.

Jacobsen, B., Bradley, C., Zidack, N., Brennemon, T., Miller, J., Washington, J., et al., 2007. Commercialization of *Bacillus mycoides* isolate BmJ as a broad spectrum biological plant disease control agent. Phytopathology 97 (7), L50.

Jain, C., Rodriguez-R, L.M., Phillippy, A.M., Konstantinidis, K.T., Aluru, S., 2018. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9 (1). https://doi.org/10.1038/s41467-018-07641-9. Article 5114.

Jallouli, W., Driss, F., Fillaveau, L., Rouis, S., 2020. Review on biopesticide production by *Bacillus thuringiensis* subsp. *kurstaki* since 1990: focus on bioprocess parameters. Process Biochem. 98, 224–232. https://doi.org/10.1016/j.procbio.2020.07.023.

Jones, D.T., Taylor, W.R., Thornton, J.M., 1992. The rapid generation of mutation data matrices from protein sequences. Bioinformatics 8 (3), 275–282. https://doi.org/ 10.1093/bioinformatics/8.3.275.

Jung, J., Jin, H., Seo, S., Jeong, M., Kim, B., Ryu, K., Oh, K., 2022. Enterotoxin genes and antibiotic susceptibility of *Bacillus cereus* isolated from garlic chives and agricultural environment. Int. J. Environ. Res. Publ. Health 19 (19). https://doi.org/10.3390/ ijerph191912159. Article 12159.

Kalyaanamoorthy, S., Minh, B.Q., Wong, T.K., Von Haeseler, A., ModelFinder, L.S. Jermiin, 2017. Fast model selection for accurate phylogenetic estimates. Nat. Methods 14 (6), 587–589. https://doi.org/10.1038/nmeth.4285.

Kardos, G., Laczkó, L., Kaszba, E., Timmer, B., Szarka, K., Prépost, E., Bányai, K., 2024. Phylogenetic analysis of the genes in D-Ala-D-Lactate synthesizing glycopeptide resistance operons: the different origins of functional and regulatory genes. Antibiotics 13 (7). https://doi.org/10.3390/antibiotics13070573. Article 573.

King, J.E., Roberts, I.S., 2016. Bacterial surfaces: front lines in host-pathogen interaction. In: Leake, M. (Ed.), *Biophysics of Infection*, Advances in Experimental Medicine and Biology, 915. Springer, Cham. https://doi.org/10.1007/978-3-319-32189-9 10.

Kuroki, R., Kawakami, K., Qin, L., Kaji, C., Watanabe, K., Kimura, Y., et al., 2009. Nosocomial bacteremia caused by biofilm-forming *Bacillus cereus* and *Bacillus*
*thuringiensis*. Intern. Med. 48 (10), 791-796. https://doi.org/10.2169/
internalmedicine.48.1885

Lazarte, J.N., Lopez, R.P., Ghiringhelli, P.D., Berón, C.M., 2018. *Bacillus wiedmannii* biovar *thuringiensis*: a specialized mosquitocidal pathogen with plasmids from diverse origins. Genome Biol. Evol. 10 (10), 2823-2833. https://doi.org/10.1093/
gbe/evy211.

Lechner, S., Mayr, R., Francis, K.P., Prüß, B.M., Kaplan, T., Wießner-Gunkel, E.L.K., et al., 1998. *Bacillus weihenstephanensis* sp. Nov. is a new psychrotolerant species of the *Bacillus cereus* group. Int. J. Syst. Evol. Microbiol. 48 (4), 1373-1382. https://doi.org/
10.1099/00207713-48-4-1373.

Leclercq, C., Allemand, P., Balcerzak, A., Branca, F., Sousa, R.F., Lartey, A., et al., 2019. A global repository for harmonized individual quantitative food consumption studies. Proc. Nutr. Soc. 78 (4), 484-495. https://doi.org/10.1017/
S0029665119000491.

Liu, Y., Lai, Q., Göker, M., Meier-Kolthoff, J.P., Wang, M., Sun, Y., et al., 2015. Genomic insights into the taxonomic status of the *Bacillus cereus* group. Sci. Rep. 5. https://
doi.org/10.1038/srep14082. Article 14082.

Liu, Y., Du, J., Lai, Q., Zeng, R., Ye, D., Xu, J., Shao, Z., 2017. Proposal of nine novel species of the *Bacillus cereus* group. Int. J. Syst. Evol. Microbiol. 67 (8), 2499-2508. https://doi.org/10.1099/ijsem.0.001821.

Ma, X., Zhang, B., Xiang, X., Li, W., Li, J., Li, Y., et al., 2024. Characterization of *Bacillus pacificus* G124 and its promoting role in plant growth and drought tolerance. Plants 13 (20). https://doi.org/10.3390/plants13202864. Article 2864.

Makuwa, S.C., Motadi, L.R., Choene, M., Liu, Y., Serepa-Dlamini, M.H., 2023. *Bacillus dicomae* sp. Nov., a new member of the *Bacillus cereus* group isolated from medicinal plant *Dicoma anomala*. Int. J. Syst. Evol. Microbiol. 73 (10). https://doi.org/
10.1099/ijsem.0.006112. Article 006112.

Mei, X., Xu, K., Yang, L., Yuan, Z., Mahillon, J., Hu, X., 2014. The genetic diversity of cereulide biosynthesis gene cluster indicates a composite transposon *Tn* ces in emetic *Bacillus weihenstephanensis*. BMC Microbiol. 14. https://doi.org/10.1186/1471-
2180-14-149. Article 149.

Messelhäusser, U., Frenzel, E., Blöchinger, C., Zucker, R., Kämpf, P., Ehling-Schulz, M., 2014. Emetic *Bacillus cereus* are more volatile than thought: recent foodborne outbreaks and prevalence studies in Bavaria (2007-2013). BioMed Res. Int. 2014.
https://doi.org/10.1155/2014/465603. Article 465603.

Miller, R.A., Beno, S.M., Kent, D.J., Carroll, L.M., Martin, N.H., Boor, K.J., Kovac, J., 2016. *Bacillus wiedmannii* sp. Nov., a psychrotolerant and cytotoxic *Bacillus cereus* group species isolated from dairy foods and dairy environments. Int. J. Syst. Evol. Microbiol. 66 (11), 4744-4753. https://doi.org/10.1099/ijsem.0.001421.

Minh, B.Q., Nguyen, M.A., Von Haeseler, A., 2013. Ultrafast approximation for phylogenetic bootstrap. Mol. Biol. Evol. 30 (5), 1188-1195. https://doi.org/
10.1093/molbev/mst024.

Minh, B.Q., Schmidt, H.A., Chernomor, O., Schremppf, D., Woodhams, M.D., Von Haeseler, A., Lan fear, R., 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37 (5), 1530-1534. https://doi.org/10.1093/molbev/msaa015. Mol. Biol. Evol., 37(5) (2020).

Mortimer, P.R., McCann, G., 1974. Food-poisoning episodes associated with *Bacillus cereus* in fried rice. Lancet 303 (7865), 1043-1045. https://doi.org/10.1016/S0140-
6736(74)90434-6.

Muigg, V., Cuénod, A., Purushothaman, S., Siegemund, M., Wittwer, M., Pflüger, V., et al., 2022. Diagnostic challenges within the *Bacillus cereus* group: finding the beast without teeth. New Microbes New Infect 49. https://doi.org/10.1016/j/
nmni.2022.101040. Article 101040.

Oren, A., Arahal, D.R., Göker, M., Moore, E.R., Rosselló-Móra, R., Sutcliffe, I.C. (Eds.), 2023. International code of nomenclature of prokaryotes. Prokaryotic code (2022 revision). Int. J. Syst. Evol. Microbiol. 73 (5a). https://doi.org/10.1099/
ijsem.0.005585. Article 005585.

Parte, A.C., Sardà Carbasse, J., Meier-Kolthoff, J.P., Reimer, L.C., Göker, M., 2020. List of prokaryotic names with Standing in Nomenclature (LPSN) moves to the DSMZ. Int. J. Syst. Evol. Microbiol. 70 (11), 5607-5612. https://doi.org/10.1099/
ijsem.0.004332.

Patino-Navarrete, R., Sanchis, V., 2017. Evolutionary processes and environmental factors underlying the genetic diversity and lifestyles of *Bacillus cereus* group bacteria. Res. Microbiol. 168 (4), 309-318. https://doi.org/10.1016/j/
resmic.2016.07.002.

Pheepakpraw, J., Kaewkod, T., Konkit, M., Krongdang, S., Jantakee, K., Praphruet, R., et al., 2023. Intraspecific diversity and pathogenicity of *Bacillus thuringiensis* isolates from an emetic illness. Toxins 15 (2). https://doi.org/10.3390/toxins15020089.
Article 89.

Priest, F.G., Barker, M., Baillie, L.W., Holmes, E.C., Maiden, M.C., 2004. Population structure and evolution of the *Bacillus cereus* group. J. Bacteriol. 186 (23), 7959-7970. https://doi.org/10.1128/jb.186.23.7959-7970.2004.

R Core Team, 2023. R: a Language and Environment for Statistical Computing R Foundation for Statistical Computing. Vienna, Austria. https://www.R-project.org/.

Rajalingam, N., Jung, J., Seo, S.M., Jin, H.S., Kim, B.E., Jeong, M.I., et al., 2022. Prevalence, distribution, enterotoxin profiles, antimicrobial resistance, and genetic diversity of *Bacillus cereus* group isolates from lettuce farms in Korea. Front. Microbiol. 13. https://doi.org/10.3389/fmicb.2022.906040. Article 906040.

Robas Mora, M., Fernández Pastrana, V.M., Oliva, L.L.G., Lobo, A.P., Jiménez Gómez, P. A., 2023. Plant growth promotion of the forage plant *Lupinus albus* Var. Orden dorado using *Pseudomonas agronomica* sp. Nov. and *Bacillus pretiosus* sp. Nov. added over a valorized agricultural biowaste. Front. Microbiol. 13. https://doi.org/
10.3389/fmicb.2022.104620. Article 1046201.

Rouzeau-Szynalski, K., Stollewerk, K., Messelhaeusser, U., Ehling-Schulz, M., 2020. Why be serious about emetic *Bacillus cereus*: cereulide production and industrial challenges. Food Microbiol. 85. https://doi.org/10.1016/j.fm.2019.103279. Article 103279.

Santé Publique France. Bulletin, 2024. Surveillance Des toxi-infections Alimentaires Collectives (TIAC). Données De la Déclaration Obligatoire. 2022. Santé publique France, pp. 1-13. https://www.santepubliquefrance.fr/maladies-et-traumatismes/
maladies-infectieuses-d-origine-alimentaire/toxi-infections-alimentaires-collecti-
ves/documents/bulletin-national/surveillance-des-toxi-infections-alimentaires-colle-
ctives-donnees-de-la-declaration-obligatoire-2022.

Sauka, D.H., Peralta, C., Pérez, M.P., Onco, M.I., Fiodor, A., Caballero, J., et al., 2022. *Bacillus toyonensis* biovar Thuringiensis: a novel entomopathogen with insecticidal activity against lepidopteran and coleopteran pests. Biol. Control 167. https://doi.org/
10.1016/j.biocontrol.2022.104838. Article 104838.

Schnepf, E., Crickmore, N., Van Rie, J., Lereclus, D., Baum, J., Feitelson, J., et al., 1998. *Bacillus thuringiensis* and its pesticidal crystal proteins. Microbiol. Mol. Biol. Rev. 62 (3), 775-806. https://doi.org/10.1128/mmbr.62.3.775-806.1998.

Schwengers, O., Jelonek, L., Dieckmann, M.A., Beyvers, S., Blom, J., Bakta, A. Goesmann, 2021. Rapid and standardized annotation of bacterial genomes via alignment-free sequence identification. Microb. Genom. 7 (11), 000685. https://doi.org/10.1099/
mgen.0.000685.

Shen, N., Yang, M., Xie, C., Pan, J., Pang, K., Zhang, H., et al., 2022. Isolation and identification of a feather degrading *Bacillus tropicus* strain Gxun-17 from marine environment and its enzyme characteristics. BMC Bio technol. 22 (1). https://doi.org/
10.1186/s12896-022-00742-w. Article 11.

Soubrier, J., Steel, M., Lee, M.S., Der Sarkissian, C., Guindon, S., Ho, S.Y., Cooper, A., 2012. The influence of rate heterogeneity among sites on the time dependence of molecular rates. Mol. Biol. Evol. 29 (11), 3345-3358. https://doi.org/10.1093/
molbev/mss140.

Stogios, P.J., Savchenko, A., 2020. Molecular mechanisms of vancomycin resistance. Protein Sci. 29 (3), 654-669. https://doi.org/10.1002/pro.3819.

Suzuki, R., Terada, Y., Shimodaira, H., 2019. Pvclust: hierarchical clustering with P-values via multiscale bootstrap resampling. http://stat.sys.i.kyoto-u.ac.jp/prog/pvclust/.

Tavaré, S., 1986. Some probabilistic and statistical problems in the analysis of DNA sequences. In: Miura, R.M. (Ed.), Some Mathematical Questions in Biology: DNA Sequence Analysis. American Mathematical Society, Providence, RI, pp. 57-86.

Thakur, P., Singh, A.K., Mukherjee, G., 2021. Isolation and characterization of alkaline pectinase productive *Bacillus tropicus* from fruit and vegetable waste dump soil. Braz. Arch. Biol. Technol. 64. https://doi.org/10.1590/1678-4324-2021200319. Article e21200319.

Tohya, M., Hishinuma, T., Watanabe, S., Shimojima, M., Ogawa, M., Tada, T., Kirikae, T., 2021. Three novel species of the *Bacillus cereus* group isolated from clinical samples in Japan. Int. J. Syst. Evol. Microbiol. 71 (9). https://doi.org/10.1099/
ijsem.0.004993. Article 004993.

Tong, C., Xiao, D., Li, Q., Gou, J., Wang, S., Zeng, Z., Xiong, W., 2024. First insights into the prevalence, genetic characteristics, and pathogenicity of *Bacillus cereus* from generations worldwide. mSphere 9 (11). https://doi.org/10.1128/msphere.00702-24. Article e00702-24.

van Tongeren, S.P., Roest, H.I., Degener, J.E., Harmsen, H.J., 2014. *Bacillus anthracis*-like bacteria and other *B. cereus* group members in a microbial community within the international space station: a challenge for rapid and easy molecular detection of virulent *B. anthracis*. PLoS One 9 (6). https://doi.org/10.1371/journal.
pone.0098871. Article e98871.

Tonkin-Hill, G., MacAlasdair, N., Ruis, C., Weimann, A., Horesh, G., Lees, J.A., et al., 2020. Producing polished prokaryotic pangenomes with the panaroo pipeline. Genome Biol. 21, 1-21. https://doi.org/10.1186/s13059-020-02090-4.

Torres Manno, M.A., Repizo, G.D., Magni, C., Dunlap, C.A., Espariz, M., 2020. The assessment of leading traits in the taxonomy of the *Bacillus cereus* group. Antonie Leeuwenhoek 113, 2223-2242. https://doi.org/10.1007/s10482-020-01494-3.

Tuipulotu, D.E., Mathur, A., Ngo, C., Man, S.M., 2021. *Bacillus cereus*: epidemiology, virulence factors, and host-pathogen interactions. Trends Microbiol. 29 (5), 458-471. https://doi.org/10.1016/j.tim.2020.09.003.

Vilain, S., Luo, Y., Hildreth, M.B., Brözel, V.S., 2006. Analysis of the life cycle of the soil saprophyte *Bacillus cereus* in liquid soil extract and in soil. Appl. Environ. Microbiol. 72 (7), 4970-4977. https://doi.org/10.1128/AEM.03076-05.

Vos, M., Didelot, X., 2008. A comparison of homologous recombination rates in bacteria and archaea. ISME J. 3 (2), 199-208. https://doi.org/10.1038/ismej.2008.93.

Wang, K., Shu, C., Soberón, M., Bravo, A., Zhang, J., 2018. Systematic characterization of *Bacillus* genetic stock center *Bacillus thuringiensis* strains using multi-locus sequence typing. J. Invertebr. Pathol. 155, 5-13. https://doi.org/10.1016/j.jip.2018.04.009.

White, H., Vos, M., Sheppard, S.K., Pascoe, B., Raymond, B., 2022. Signatures of selection in core and accessory genomes indicate different ecological drivers of diversification among *Bacillus cereus* clades. Mol. Ecol. 31 (13), 3584-3597. https://doi.org/
10.1111/mec.16490.

Wickham, H., 2016. ggplot2: Elegant Graphics for Data Analysis. Springer, New York, NY.

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al., 2021. ClusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation 2 (3). https://doi.org/10.1016/j.xinn.2021.100141.

Xu, S., Dai, Z., Guo, P., Fu, X., Liu, S., Zhou, L., et al., 2021. Compact visualization of richly annotated phylogenetic data. Mol. Biol. Evol. 38 (9), 4039-4042. https://doi.org/
10.1093/molbev/mnsab166.
Yang, Z., 1995. A space-time process model for the evolution of DNA sequences. *Genetics* **139**(2), 993–1005. https://doi.org/10.1093/genetics/139.2.993.

Yu, G., 2022. Data Integration, Manipulation and Visualization of Phylogenetic Trees.
Chapman and Hall/CRC. ISBN: 9781032233574.

Yu, G., Smith, D.K., Zhu, H., Guan, Y., Lam, T.T.Y., 2017. An R package for visualization and annotation of phylogenetic trees with their covariates and other associated data.

Methods Ecol. Evol. 8 (1), 28–36. https://doi.org/10.1111/2041-210X.12628.