A total of 100 samples of tomatoes, bell peppers, and salads from a central retail distribution facility in Switzerland were screened for B. cereus sensu lato (Table 1 and Supplementary Table 1). 10 g of each sample was added to 90 g peptone salt solution (0.1% enzymatic casein digest) and homogenized in a stomacher@ 400 circulator (Seward) for 60 s. Subsequently, B. cereus s.l. counts were determined by plating in serial dilutions on selective media (Mannitol-Egg Yolk-Polymyxin [MYP] agar) and incubated at 30°C overnight. Up to five presumptive B. cereus s.l. colonies (mannitol negative, phospholipase C positive phenotype) per sample were selected and whenever different colony morphologies were observed, these were included in the selection. Isolates were purified and preserved in glycerol stocks at −80°C until further characterization.

Pure cultures of all isolates were grown on T3 agar (Travers et al., 1987) for 3 days at 30°C. A tiny amount of colony material was resuspended in 10 μL sterile water and screened for crystals of diamond, spherical or bipyramidal shape using phase-contrast microscopy (EFSA BIOHAZ Panel, 2016). Isolates exhibiting characteristic parasporal crystals were defined as B. thuringiensis.

All B. thuringiensis were screened by FT-IR spectroscopy to detect clusters of related isolates and select a subset of strains for whole-genome sequencing. To this end, metabolic fingerprints of all isolates were generated as previously described (Ehling-Schulz et al., 2005). Briefly, isolates were cultivated as lawns on tryptone soy agar (TSA) (Oxoid) at 25°C for 24 h. Subsequently, one loop of cell material was suspended in 300 μL sterile deionized water, subjected to ultrasonication for 5 × 1 s at 100% power with a Bandelin Sonopuls HD2200 (Bandelin electronic), and the bacterial suspension was spotted on a zinc selenite (ZnSe) optical plate. After drying the plates, infrared absorption spectra were recorded using an HTS-XT microplate adapter coupled to a Tensor 27 FT-IR spectrometer (Bruker Optics). Spectra in the range of 4,000–500 cm^–1 were acquired in transmission mode with the following parameters: 6 cm^–1 spectral resolution, zero-filling factor 4, Blackmann-Harris 3-term apodization, and 32 interferograms were averaged with background subtraction for each spectrum. The quality of FTIR spectral data was evaluated using the OPUS software (version 7.5; Bruker Optics). Second derivatives were calculated using the Savitzky-Golay algorithm with 11 smoothing points and the derivative spectra were unit vector normalized for further data processing. The so-called fingerprint region (spectral region of 1,500–800 cm^–1) was chosen for chemometric analyses. Hierarchical cluster analysis (HCA) was carried out using Ward’s algorithm.

For short-read sequencing, genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen). Libraries were prepared using the Nextera DNA Flex Library Preparation Kit (Illumina), and sequencing was performed on the Illumina MiniSeq platform with 2 × 150 bp paired-end chemistries. Illumina read adapters and low-quality bases were trimmed with TrimGalore v0.6.6¹ and quality metrics computed using FastQC v0.11.9². Draft genomes were assembled using SPAdes v3.14.1 (Bankevich et al., 2012) implemented in shovill v1.1.0³. Additionally, all seven biopesticide isolates obtained during this study were long-read sequenced. For nanopore sequencing, genomic DNA was extracted using the MasterPure Complete DNA and RNA Purification Kit (Lucigen). Multiplex libraries were prepared using the SQK-LSK109 ligation sequencing kit with the EXP-NBD114 native barcoding expansion kit (Oxford Nanopore Technologies). Libraries were sequenced on a MinION Mk1B device using the FLO-MIN106 (R9) flow cell (Oxford Nanopore Technologies). Adapters were trimmed using Porechop v0.2.4⁴ and quality assessed with LongQC v1.2.0 (Fukasawa et al., 2020). Hybrid assemblies were produced with Unicycler 0.4.8 (Wick et al., 2017). Assembly quality metrics were assessed using QUAST v5.0.3 (Gurevich et al., 2013).

Multi-locus sequence types (MLST), panC types, and the presence of virulence genes were determined using BTyper3 and srst2 0.2.0 (Inouye et al., 2014; Carroll et al., 2020a). The completeness of virulence genes and their promoter regions was additionally investigated using BLAST (Altschul et al., 1990) implemented in abricate 1.0.1⁵ and putative disruptions manually inspected in assemblies annotated with prokka v1.14.6 (Seemann, 2014). Prophages were detected using PHASTER (Arndt et al., 2016). sph genes disrupted by the insertion of a BceA1-like bacteriophage were identified by screening genomes for the sequence of open reading frame BG08_RS31255 (strain HD-1, GCF_000835235.1;>99% alignment coverage,>99% sequence identity in BLAST), which contains elements of sph and of a site-specific integrase from the bacteriophage.

Core genome SNPs (cgSNPs) among CC8 and CC23 isolates were detected from Illumina read data using the CFSAN SNP pipeline v2.2.1 (Davis et al., 2015) with chromosomes of strains HD-1 (CC8, GCF_000835235.1), and BGSC 4Q7rifR (CC23, GCF_013113775.1) as references, respectively. Phages, IS elements, and repeat regions in reference chromosomes were identified using phastaf v0.1.0⁶, ISEScan v1.7.2.3 (Xie and Tang, 2017), and NUCmer v3.1 (Marçais et al., 2018) and masked prior to read-mapping. cgSNP, i.e., SNP sites present in all genomes, were extracted from the SNP alignment by removing missing sites using SNP-sites v2.5.1 (Page et al., 2016), and core-genome SNP (cgSNP) distances were determined with snp-dists v0.7.0⁷. Maximum-likelihood phylogenetic trees were constructed from the SNP alignment using IQ-TREE v2.0.3 (Nguyen et al., 2015) with the generalized time-reversible (GTR) model and gamma distribution with 100 bootstraps. The number of invariant sites was estimated from core genome alignments generated with parsnp v1.5.3 (Treangen et al., 2014) and passed to IQ-TREE. A phylogenetic tree of all 33 isolates was created using IQ-TREE directly from a core genome alignment generated with parsnp.

Whole-genome SNPs (wgSNPs) were detected by mapping reads of food/human isolates to the assembly of the phylogenetically closest biopesticide isolate (without any masking) using the CFSAN SNP pipeline. Read data of the respective reference biopesticide isolate were included in the pipeline as a reference to calculate wgSNP distances. The percentage of reference assembly bases covered by sequencing reads was determined with BBMap⁸.

A total of 100 food samples (tomatoes, bell pepper, endives, rocket salad, baby spinach, and baby leaf salad) were collected at the retail level. Presumptive B. cereus s. l. isolates were identified in 27 samples. From these 27 samples, a total of 78 isolates (up to 5 isolates per sample) were screened by phase-contrast microscopy. A crystal-producing phenotype, characteristic for B. thuringiensis, was identified for 49 isolates (63%) from 14 food samples (14%). Positive samples included tomatoes (7/39, 18%), bell pepper (4/21, 19%), rocket salad (2/11, 18%), and endives (1/18, 6%) (Table 1). No B. thuringiensis was detected in samples from baby spinach (n = 7) and baby leaf salad (n = 4). Remarkably, all investigated B. cereus s. l. isolates from tomatoes (n = 26), bell pepper (n = 12), and rocket salad (n = 8) were identified as B. thuringiensis. In the 14 samples positive for B. thuringiensis, presumptive B. cereus counts ranged from 100 to 40,000 CFU/g and were thus below the alert threshold of 10⁵ CFU/g (Table 1 and Supplementary Table 1). In 13 food samples, only non-crystal-producing B. cereus s. l. isolates were detected. These were particularly common in endives (9/18, 50%) and baby leaf salad (2/4, 50%) and reached bacterial loads of up to 1,800 CFU/g.

B. cereus s. l. isolates from all food samples were subjected to FT-IR spectroscopy for dereplication and selection of representative isolates to be included in the WGS analysis. B. thuringiensis isolates from the biopesticide products Agree, B401, Delfin, DiPel, Novodor, Solbac, and XenTari were included as references. Closely related isolates were grouped into clusters using Ward’s hierarchical clustering method. Out of each cluster containing B. thuringiensis, either all (for small clusters and singletons) or a minimum of three isolates from distinct samples were subjected to WGS analysis for further characterization. Two samples were represented by two isolates each which fell into different FT-IR clusters (Supplementary Table 2).

Overall, whole-genome sequences of 13 biopesticide product isolates, 15 food isolates, and 5 isolates linked to two foodborne outbreaks were investigated to identify potential matches of food/outbreak isolates with biopesticide strains. Of these 33 B. thuringiensis isolates, 27 were sequenced as part of this study (Table 2). These included 9 food isolates obtained from 7 samples as part of this study. To obtain a more comprehensive picture of possible links between biopesticide strains, food isolates, and outbreaks, the collection was supplemented with 18 previously described isolates originating from biopesticide products (n = 7), food (n = 9), or human fecal samples (n = 2) (Johler et al., 2018). Three of these food isolates and the two fecal isolates were linked to foodborne outbreaks of diarrheal disease (Johler et al., 2018); three food isolates were recovered during controls by a cantonal laboratory (CH_69, CH_72, CH_81) from lasagna, vegetable juice, and sauce; two isolates were obtained from honey during the production process; and one isolate was obtained from tarragon at the retail level. Furthermore, publicly available WGS data of 6 previously sequenced isolates from additional biopesticide products (Bonis et al., 2021) were included for genomic analyses.

All 33 genomes contained the B. thuringiensis-defining cry genes and were assigned to panC group IV. The 5 outbreak-associated isolates and 15 food isolates belonged to one of two closely related lineages ST8 (subspecies kurstaki) or ST15 (subspecies aizawai), which together form CC8. Most (11 out of 15) biopesticide isolates investigated in our study also fell into these two lineages. SNP analyses suggested that all food/fecal isolates originated from one of six distinct biopesticide strains: 19 of the 20 isolates differed by ≤ 4 wgSNPs (≤ 2 cgSNPs) from their phylogenetically closest biopesticide isolate (Figure 1) and food isolate CH_65 differed by 6 wgSNPs (2 cgSNPs) from its closest biopesticide isolate. The food isolates matched with biopesticide strains ABTS-351, ABTS-1857, ABTS-1857/B401, GC-91, SA-11, or SA-12. The five outbreak-associated isolates clustered in a sublineage of ST15 and differed by 0–3 wgSNPs (0–2 cgSNPs) from the biopesticide strain ABTS-1857, suggesting this biopesticide as a highly likely origin. Whereas most biopesticides are applied on a wide spectrum of crops, the biopesticide B401 is specifically recommended for the control of wax moths on honeycombs. Accordingly, the two isolates recovered from honey products corresponded to the B401 isolate (0 and 2 wgSNP/cgSNP difference, respectively).

The toxin gene profile varied between biopesticide isolates (Table 3). All biopesticide isolates except CH_187 (NB-176) carried intact genes and associated 5’ intergenic regions encoding cytotoxin CytK2 (cytK2), hemolysin BL (hblA, hblB, hblC, hblD), and non-hemolytic enterotoxin (nheA, nheB, nheC). In CH_187 (NB-176), cytK was absent and the nheA promoter disrupted by a transposase inserted 13 bp upstream of the translation initiation site, likely preventing effective transcription. As expected for B. thuringiensis, none of the isolates possessed the ces locus (encoding the emetic toxin cereulide).

In multiple B. thuringiensis genomes, the sph gene (sphingomyelinase) was disrupted by the insertion of prophages (Figure 2 and Table 3). Hybrid assemblies from biopesticide isolates CH_183 (ABTS-351, ST8) and CH_164 (SA-11, ST8) revealed the presence of the intact 42 kb Bacillus prophage BtCS33 (NC_018085.1, >99% alignment coverage and identity) inserted into the 3′-terminus of sph. Long-read assemblies from isolate CH_186 (GC-91, ST15) generated using Canu or Raven contained the intact 42 kb Bacillus prophage BceA1 (GCF_003047795.1, > 99% alignment coverage and identity) inserted into the 3′-terminus of sph. By contrast, sph was located at the end of a split contig in a hybrid assembly of CH_186 generated with Unicycler. Verification by PCR suggested that both intact and phage-disrupted sph loci are detectable in CH_186, possibly resulting from subpopulations of lysogenic, pseudo-lysogenic, lytic, and/or naïve cells. In sequencing read mapping analyses, the average base coverage in prophage regions of CH_186 was > 4-fold higher than in flanking regions, further indicating the presence of extrachromosomal DNA of active bacteriophages. Almost half of the BtCS33 and BceA1 genomes are homologous, including the identical site-specific integrases GP25 and P30, respectively (Figure 2). A 19 bp identical region in sph and downstream of the integrase genes was identified as the insertion site. In short-read assemblies, phage insertions could be identified by an incomplete sph sequence (96.85% alignment coverage). Accordingly, disrupted sph was also detected in all other biopesticide isolates from ST8, whereas biopesticide isolates from CC23 and two out of three biopesticide isolates from ST15 harbored intact sph genes. The toxin gene profile (sph completeness and nhe, hbl, cytK, ces) of food and fecal isolates always corresponded to their respective closest biopesticide isolate (shown in Figure 1).