Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France

Bacillus thuringiensis (Bt) belongs to the Bacillus cereus (Bc) group, well known as an etiological agent of foodborne outbreaks (FBOs). Bt distinguishes itself from other Bc by its ability to synthesize insecticidal crystals. However, the search for these crystals is not routinely performed in food safety or clinical investigation, and the actual involvement of Bt in the occurrence of FBOs is not known. In the present study, we reveal that Bt was detected in the context of 49 FBOs declared in France between 2007 and 2017. In 19 of these FBOs, Bt was the only microorganism detected, making it the most likely causal agent. Searching for its putative origin of contamination, we noticed that more than 50% of Bt isolates were collected from dishes containing raw vegetables, in particular tomatoes (48%). Moreover, the genomic characterization of isolates showed that most FBO-associated Bt isolates exhibited a quantified genomic proximity to Bt strains, used as biopesticides, especially those from subspecies aizawai and kurstaki. Taken together, these results strengthen the hypothesis of an agricultural origin for the Bt contamination and call for further investigations on Bt pesticides.


Introduction
Bacillus thuringiensis (Bt) is a ubiquitous, Gram-positive bacterium, well known for its ability to produce insecticidal crystals during sporulation [1,2]. These crystals are composed of δendotoxins classified within a large family of Cry/Cyt toxins (more than 850 halotypes described so far), firstly according to their specificity of action against insects [3], and secondly according to their sequence similarity [4]. When activated in the midgut of larvae, the toxins The aim of the present study was first to estimate the frequency of Bt isolation in the context of FBOs that occurred in France and second, to assess their genetic diversity and to compare them with Bt strains isolated from commercially available products. Bt isolates were identified from a large collection of FBO-associated Bc isolates from France, based on the detection of parasporal crystals. Phenotypic and genomic characterizations were then conducted to compare isolates both pairwise and against a representative panel of Bt strains from the European market.

FBO-associated isolates and data collection
For this study, we used a total of 1,263 B. cereus sensu lato (Bc) isolates previously collected in France from food during investigations of 250 French foodborne outbreaks that occurred in the country between 2007 and 2017. Outbreaks were declared to local health authorities, who collected epidemiological data through interviews or questionnaires. They included the date of the FBO, the nature of the consumed dishes, the number of patients, the symptoms, and the time before onset of symptoms. Based on these data, hypotheses about the etiological agent were established and food samples were collected for analysis. The detection and the counting of Bc isolates from food were performed by certified French laboratories, as previously described [29], and according to the International Organization for Standardization (ISO) 7932 standard method. Between 1 and 15 colonies per FBO or foodstuff were selected and transferred to the "Staphylococcus, Bacillus and Clostrium" unit of ANSES for further analysis and characterization. Information related to the putative detection of other food microorganisms was also collected.

Identification of Bt isolates within the Bc collection
The FBO-associated Bc collection was screened for the presence of Bt isolates. As no standardized procedure was available to distinguish Bt from other Bc, we used a method consensually accepted at the international standardization level (ISO/TC34/SC9/WG20), based on the microscopic detection of parasporal insecticide crystals, as follows. Bc isolates were cultivated overnight on tryptone soya agar yeast extract (TSAYE) agar plates (Biokar) at 30˚C, then on the solid sporulating medium hydrolysate of casein tryptone (HCT), supplemented with 0.3% glucose, at 30˚C [30]. A sample of the culture was then resuspended in 10 μL of distilled water and screened by phase-contrast microscopy every 24h until 72h of incubation. Images were acquired with an Olympus BX51, using the Archimed software. When crystals were observed, corresponding Bc isolates were assigned to the Bt species.

Examination of incriminated foodstuffs
To investigate the origin of contamination for the Bt isolates, we examined the nature of the food products from which they were isolated, comparatively to Bc (non-Bt) isolates. Since the Bc collection exhibited a clonal redundancy per FBO, 59 Bt and 437 non-Bt representative isolates were selected, based on their original phenotypic and genotypic profiles (see respectively parts "Phenotypic characterization of Bt isolates" and Genotypic characterization of Bt isolates) per FBO and per foodstuff. Then we calculated the frequency of Bt/non-Bt isolates association with eight types of foodstuff: tomatoes (including cooked tomatoes but excluding tomato sauce), lettuce, raw vegetables, fruits (raw and cooked), starch products, meat, fish or seafood and spices or dried herbs. A panel of 19 phytopharmaceutical and biocides products containing single, non-genetically  modified Bt strains was collected (Table 2). These preparations were selected to cover all the biopesticide uses of Bt authorized in France under Regulations (EC) No 1107/2009 and No 528/2012. In agreement with the labeling satisfying the requirements laid down in the Regulation (EC) No 1107/2009, they included altogether 10 different Bt strains, divided into 4 subspecies: aizawai, kurstaki, israelensis and morrisoni. For isolation, one g of each product was resuspended in 9 mL of triptic soy broth (TS) until complete dissolution, and bacteria were isolated by successive plating on selective media (Mossel, Biomérieux).

Phenotypic characterization of Bt isolates
The lecithinase and hemolytic activities of FBO and insecticide Bt isolates were assessed according to the ISO 7932 standard method, as follows. The presence or absence of lecithinase activity was detected by bacterial cultivation on Mossel agar plates (Biomérieux) for 18 to 24 hours (h) at 30˚C. When expressing the lecithinase enzyme, isolates exhibited an opaque halo of precipitation surrounding the colonies due to the presence of an egg yolk emulsion. Hemolytic activity was evaluated after inoculating Columbia medium agar plates (Biomérieux) containing 5% sheep blood, with single colonies from TSAYE agar plates 24h-cultures. Plates were incubated for 18 to 24 hours at 30˚C before assessing the presence of the hemolytic halo.

DNA isolation for genotypic characterization
The DNA of Bt isolates was extracted after overnight cultures at 30˚C on TSAYE agar plates, using a DNeasy Blood and Tissue Kit (Qiagen). Briefly, a few colonies were resuspended in 180 μL of 1M Tris-HCl, 0.25M Sodium-EDTA, 1.2% TritonX-100 and 2% lysozyme, and incubated for 1 hour at 37˚C. The following steps comply with the manufacturer's recommendations. DNA was quantified by absorbance measurement at 260 nm, using the NanoDrop1000 spectrophotometer (Thermo Fisher Scientific).

Genotypic characterization of Bt isolates
Bt isolates were genotypically characterized according to the method detailed in S1 Table. The presence of the cytK1/2, ces, hlyII, nheA/B/C and hblA/C/D genes was determined by PCR with the ProFlex PCR System (Applied Biosystems), using respective primers and methods (S1 Table). To assign Bt isolates to one of the seven described phylogenetic groups, a part of the panC gene was amplified as described (S1 Table) and sequenced (Eurofins MWG Operon). For publicly available sequences, the full sequence of panC was extracted with BLAST [31], using the ATCC 14579 panC sequence (AE016877.1) as the query. The phylogenetic classification was then performed online (https://www.tools.symprevius.org/Bcereus/) using a public algorithm based on panC sequence similarity and the statistical significance of matches with database sequences [9]. All the PCR products were analyzed by electrophoresis onto 2% agarose gels (made with a 50/50 mix of Seakem1 GTG™ Agarose and NuSiev1 GTG™ Agarose, Lonza), with a run in 1X TBE, for 1.5 hours at 90V.

M13-PCR typing
To study the diversity of Bc isolates, we used a coliphage M13 sequence-based PCR (M13-PCR) derived from a Random Amplified Polymorphic DNA (RAPD) technique and adapted from [32], according to the method described in S1 Table. The PCR products were separated by electrophoresis onto 1% agarose gels (made with a 50/50 mix of Seakem1 GTG™ Agarose and NuSiev1 GTG™ Agarose, Lonza), with a run in 1X TBE, for 10 min at 50V, followed by 3.5h at 90V. Gels were stained with ethidium bromide. For this study, the electrophoretic patterns of 59 representative FBO-Bt and 19 representative insecticide Bt isolates were analyzed and compared using Bionumerics 7.6. A dendrogram was constructed based on pairwise Dice similarity coefficients calculations [33] and UPGMA clustering, with tolerance and optimisation set at 1%. The Bt isolates were tested three times by M13-PCR typing and similar visual patterns were observed.

Detection of the Nhe and Hbl enterotoxins
To confirm the in vitro expression of nhe and hbl, the toxins production was assessed from a panel of 10 representative commercial Bt isolates and 21 representative Bt isolates, selected for their original phenotypic and genotypic profiles per foodstuff and per FBO from group A (i.e. for which no other microorganism could be detected). Each Bt isolate was cultivated on TSAYE agar plates for 18 hours at 30˚C, then in 10 mL of brain heart infusion broth (BHI, Biomérieux) inoculated with a single colony, for 6 hours at 37˚C with stirring. Production of the Nhe and Hbl enterotoxins was then estimated with two immunological tests, using the BDE VIATM kit (3M-Tecra) and the BCET-RPLA kit (Oxoïd), respectively, as previously described [34]. The value indices obtained for FBO-Bt and commercial Bt were compared using the nonparametric Wilcoxon statistical test.

Whole genome sequencing (WGS)
The DNA of Bt/Bc isolates was extracted using the KingFisher Cell and Tissue DNA kit (Ther-moFisher), as follows. Bacteria were grown overnight at 30˚C, successively on TSAYE agar plates and then in 10 mL of BHI. Cells from 2 mL aliquots were lysed in 600 μL of 40 mM EDTA, 2 mg/mL lysozyme for 1 hour at 37˚C, before being resuspended in 200 μL of lysis buffer and 25 μL of proteinase K, and incubated for 15 minutes at 70˚C then 30 minutes at 80˚C. RNA was removed by incubating for 5 minutes at 37˚C with 12 μg/mL of RNAse A (Promega). The following purification steps were performed using the KingFisher Duo Prime system, according to the manufacturer's recommendations. Final elution was performed in 100 μL of PCR grade water. To estimate the DNA purity, we measured the absorbance of DNA at 230, 260 and 280 nm using the NanoDrop2000 (ThermoScientific). The ratio A260/A280 and A260/A230 ratios were expected to be comprised between 1.5 and 2.5. The concentration of the solutions was measured with the Qubit Fluorometer (Invitrogen) and the dsDNA HS Assay kit (Invitrogen), according to the manufacturer's recommendations. The integrity of genomic DNA was checked by visualization onto 0.8% agarose gels (made with a 50/50 mix of Seakem1 GTG™ Agarose and NuSiev1 GTG™ Agarose, Lonza), with a run in 1X TBE, for approximately 2 hours at 90V. Isolated DNA was sent to the ICM (« Institut du Cerveau et de la Moëlle épinière») for paired-end read sequencing (2x150bp), using the Nextera XT DNA Library Prep Kit and the Nextseq500 sequencing system (Illumina).

Genome assembly
The reads were assembled with the in-house workflow ARtWORK v1.0 [35], with default parameters. After read decompressing, short paired-end reads were mapped against the closely related reference genome identified by estimating the Jaccard index with Mash v2.0 [36], among a selected collection of high-quality fully closed assemblies, representative of the species (S4 Table). The samples presenting a depth of coverage against the reference lower than 50X were excluded, those between 50X and 100X were retained, and those higher than 100X were normalized at 100X with Bbnorm v38.22 (https://jgi.doe.gov/data-and-tools/bbtools/). Reads were trimmed with Trimmomatic v0.38 (i.e. adapters sequences, quality phred score of 30 and minimum length of 50 bp) and assembled using SPAdes v3.13.0, with the careful option to reduce the number of mismatches and short indels and the automatic k-mer sizes selection [37]. The closest reference genome was used to perform reference-based scaffolding with MeDuSa v1.6 [18,38] and gap filling with GMcloser v1.5 [39]. Scaffolds smaller than 150 bp were removed from the final assembly.

k-mer phylogeny
A phylogenetic tree of 234 Bc/Bt genome assemblies (listed in S4 Table) was built from the number of shared k-mers (kmer size of 17 and sketch size of 10,000) with the in-house script QuickPhylo v1.0 (https://github.com/afelten-Anses/QuickPhylo). This tool used Mash to compute a distance matrix based on the Jaccard index and cluster sequences with the neighborjoining (NJ) method. The midpoint approach was used to root the tree and its visualization was obtained with iTOL v5 [40]. Sixty-two publicly available genome sequences were added from NCBI to enhance the genetic diversity of the dataset.

Single nucleotide polymorphism (SNP) analysis
To estimate the genetic distance separating FBO and insecticide isolates, we conducted an SNP analysis using the iVARCall2 v1.0 workflow [41] from sequence reads of 154 Bc isolates belonging to phylogenetic group IV (listed in S4 Table). These samples comprised 19 insecticide Bt isolates and 135 FBO isolates (123 Bt and 12 non-Bt). The closed genome of Bt strain HD73 was used as a reference for read mapping (accession NC_020238.1). The IVARCall2 workflow generated a matrix of pairwise SNP distances for each pair and pseudogenome sequences. These pseudogenomes were reconstructed by replacing the genotypes of detected variants into the reference genome for each paired-end reads of the dataset. Recombination events were excluded to infer the presented phylogenomic reconstruction using ClonalFra-meML v1.12, run with 100 simulations [42]. A phylogenetic reconstruction was performed from pseudogenomes using IQ-TREE v1.6.9 [43]. The K3P model was selected as the best fitting model of reconstruction and 10,000 ultrafast bootstraps were generated. Visualization of trees was obtained with iTOL v5 [40].

Comparison of pairwise SNP distances
The comparison of pairwise SNP distances was performed from the matrix generated by iVARCall2 [41], using an R script adapted from a public script (https://github.com/andersgs/ harrietr). According to the phylogenetic reconstruction and calculated SNP distances, each isolate was attributed to one of the 5 defined clusters (hereafter referenced from "a" to "e"). Then, intra-and inter-clusters SNP distances were computed and compared with a Wilcoxon statistical test. Their distribution was illustrated with box plot diagrams.

Accessory genome analysis
As an addition to the SNP analysis approach, the accessory genomes of the 154 FBO and insecticide isolates were compared. A tree was built using ROARY v3.12.0 with default parameters [44] and iTOL v5 [40]. Similarly to the comparison of pairwise SNP distances, the intra-and inter-clusters pairwise numbers of different genes were computed and compared.

Distinction of Bacillus thuringiensis (Bt) among FBO-associated Bc isolates
A collection of 1,263 Bc isolates associated with 250 distinct FBOs was established and characterized between 2007 and 2017. Retrospectively, we tested the ability of these isolates to produce crystals under sporulating conditions by phase-contrast microscopy, in order to identify the putative presence of Bt isolates. Similarly to the positive control CIP53137, Bt isolates exhibited the production of intracellular inclusions after 48 h of incubation, followed by the release of crystals and spores (S1 Fig). By contrast, the Bc strain ATCC14579, used as a negative control, was unable to produce crystals. Hence, we identified 143 Bt isolates, representing 11.3% of the total collection of 1,263 Bc. More importantly, these 143 Bt isolates were associated with 49 toxic episodes (listed in Table 1), representing 19.6% of the 250 FBOs considered in the present study.

Other food pathogens in Bt-associated FBOs
To determine whether Bt isolates were the putative etiological agents of the 49 FBOs, we examined the potential co-occurrence of other putative foodborne pathogens classically screened for during FBO investigations: Salmonella spp, Staphylococcus aureus, Clostridium perfringens and enteric viruses when symptoms occurred more than 24 h after ingestion of suspected food. The collected data revealed that for 19/49 FBOs (38.8%), no other food pathogen could be detected (Table 1). This ranks Bt isolates first among the hypothetical etiological agents of these 19 FBOs, on the basis of current knowledge. These FBOs were attributed to group A for this study, the others being attributed to group B.

Epidemiology of Bt-associated FBOs
According to the available epidemiologic data, more than 673 patients were affected by the 49 Bt-associated FBOs, and more than 330 patients when restricted to the 19 Bt-FBOs from group A. The symptoms reported were of the gastroenteritis type, such as diarrhea, vomiting, nausea, and abdominal pain, with a median onset time of 5 h (from 15 minutes to more than 24 h). Depending of the FBO, the rate of Bt contamination ranged from 1x10 2 to more than 1x10 7 CFU/g of food (median of 9x10 2 ). In 44/49 FBOs (89.8%), the level of Bt contamination was found to be below the alert threshold of 1x10 5 CFU/g established by the French Ministry of Agriculture [16].

Bt versus non-Bt isolates in foodstuffs
To search for food products specifically associated with Bt isolates, we compared the frequency of isolation of representative FBO-Bt and -Bc (non-Bt) isolates, from dishes containing eight types of food products (Fig 1). Our data revealed that Bt was isolated more frequently than non-Bt in the presence of raw vegetables, in particular tomatoes and lettuce (respectively 7, 9 and 15 times more), and in the presence of fruits (3 times more). By contrast, Bt was found to be less commonly associated than other Bc to starch products, meat, fish or seafood, and spices or dried herbs.

FBO-associated Bt isolates and pesticide Bt strains
To explore the putative path of agricultural Bt contamination via pesticides, a panel of 19 phytopharmaceutical products formulated with 10 different Bt strains was collected ( Table 2). After isolating their bacterial content on specific media, we compared them to FBO-associated Bt, from a selection of representative isolates (i.e. exhibiting original phenotypic and genotypic profiles per FBO and per foodstuff ( Table 3, S2 Table). Almost all Bt isolates were found to be hemolytic and exhibited a lecithinase activity on Mossel plates, with the exception of one FBO-Bt (14SBCL08) and one commercial Bt (18SBCL485), found to be lecithinase-negative. All commercial Bt belonged to phylogenetic group IV, as well as 98% of FBO-Bt. We then searched, by conventional PCR, for the presence of genes encoding six virulence factors. The majority of Bt isolates (94%) exhibited a common signature, defined by the presence of cytK2, nheA/B/C and hblC/D/A, and the absence of ces, hlyII and cytK1. By contrast, the commercial strain NB-176 (isolate 18SBCL485) was cytK2-negative, and all commercial Bt ssp. israelensis and morrisoni strains were found to be hlyII-positive, as well as two FBO-Bt isolates (08CEB128BAC and 17SBCL429).

PLOS ONE
The lecithinase and hemolysis activities were searched on specific media as described in materials and methods. The presence of virulence genes cytK1, cytK2, ces, hlyII, nhe (A/B/C) and hbl (C/D/A) was determined by conventional PCR. The assignment to phylogenetic group IV was realized by partial sequencing of the panC gene. The percentages were calculated from panels of 59 representative FBO-Bt isolates and 19 insecticide strains. The detailed characterization of all Bt isolates is listed in S2 Table.

In vitro production of Nhe and Hbl
After detecting the presence of nhe and hbl for most of Bt isolates, we assessed their expression in vitro by detection of toxin components in culture supernatants. The experiment was performed with a panel of 31 isolates (composed of 21 Bt isolated from 19 FBOs from group A and 10 commercial Bt), using two immunological tests (Fig 2, S3 Table). Results indicated that the Nhe production indices of Bt isolates from FBO and pesticides were analogous (i.e. not significantly different, p-value = 0.9003), with respective Nhe absorbance medians of 0.425 and 0.405. Likewise, the Hbl production levels of FBO-Bt isolates was found to be comparable to the one of pesticide Bt (median Hbl dilution = 1/64 in both cases, p-value = 0.4168). The  18SBCL483A isolate (commercial Bt strain BMP144) was the only Bt isolate found to be associated with a level of Hbl production below the detection limit in tested conditions.

M13-PCR typing
With the aim of comparing Bt isolates more precisely, they were analyzed by M13-PCR typing and Dice's coefficient similarity calculations (Fig 3). Based on a similarity threshold of 80%, commonly used for this type of method with Bc [32,45], Bt isolates were clustered into three groups named 1 to 3. Regarding the commercial strains, Bt isolates were grouped according to their subspecies, as follows. Group

Genomic relationships between Bt isolates from FBO and commercial products
To further compare Bt isolates, we carried out analyses of their whole genomes. The DNA of 172 isolates (i.e. 152 FBO-associated strains, 19 insecticide strains and one reference strain) was purified and sequenced before genome assembly. To evaluate their genomic relatedness within the Bc group, a phylogenetic tree based on shared k-mers was built, together with 62 publicly available sequences of various Bc strains (Fig 4). Most genomes were found to be organized according to their phylogenetic group based on panC sequencing; groups I and IV were the most distant., Interestingly, the k-mer phylogeny revealed a dominant clade (highlighted in grey) composed of 144 closely related genomes of group IV Bt isolates. Except for two strains (08CEB128 and 17SBCL429), this clade included all the FBO-associated Bt analyzed in the present study. They included all the collected commercial strains from ssp. aizawai and kurstaki. By contrast, commercial Bt isolates from ssp. israelensis and morrisoni were located out of this clade.

Main clusters of FBO and commercial Bt isolates
To refine and consolidate the k-mer-based phylogeny, we performed an SNP-based analysis of the core genome of 154 isolates from group IV, comprising all commercial Bt isolates and all FBO-associated Bt isolates (except for 08CEB128 because of its large genetic distance). From a new phylogenetic reconstruction (Fig 5A) and the computed pairwise SNP distances (not shown), we highlighted the presence of 5 clusters (named a to e). Cluster e, that grouped together 17 isolates distant from the rest of the dataset, was removed from the tree for better visualization (Fig 5B). In agreement with the k-mer phylogeny, commercial Bt isolates aizawai and kurstaki were scattered among Bt isolated from 47 toxic episodes, and distributed into clusters a to d. Clusters a and b comprised isolates related to 9 and 11 FBOs, respectively along with commercial Bt strains GC-91 (cluster a), and ABTS-1857 (cluster b). Cluster c included 13 FBOs with pesticides SA-12, PB-54 and EG2348. Lastly, the commercial strains SA-11 and ABTS-351 were located in cluster d, as well as Bt associated with 21 FBOs (Fig 6). In the case of 4 outbreaks (23, 33, 36 and 45), Bt isolates were dispatched into 2 distinct clusters, suggesting

Fig 3. Dendrogram of representative FBO-and commercial Bt isolates based on M13-PCR typing data.
After migration of M13-PCR products onto 1% agarose gels, the acquired images were analyzed with Bionumerics 7.6. The dendrogram was constructed using Dice coefficients and UPGMA clustering, with tolerance and optimisation set at 1%. Three groups of isolates were defined (1 to 3) based on a calculated similarity level above 80%.
https://doi.org/10.1371/journal.pone.0246885.g003 from pesticides, and 32 from the environment (soils or insects). Genomic sequences were obtained in this study after genomic DNA extraction and NGS for 172 isolates, or collected from the NCBI database for 62 isolates. Detailed information related to each isolate is listed in S4 Table. The phylogenic tree was obtained with an in-house script using Mash as distance estimation and neighbor-joining (NJ) as the clustering method. Tree visualization was performed with iTOL v5 [46]. The IDs of Bt strains isolated from biopesticides were highlighted according to their respective subspecies (green: morrisoni, blue: israelensis, yellow: kurstaki, and orange: aizawai). The genome IDs from the dominant clade was highlighted in grey. https://doi.org/10.1371/journal.pone.0246885.g004

PLOS ONE
the presence of two Bt populations. Moreover, this distribution was confirmed with a phylogenetic reconstruction built from clustering of the accessory genomes (S2 Fig). The pairwise SNP distances between Bt isolates ranged from 0 to 10 SNP within a same cluster, with median distances of 2, 0, 1 and 1 SNP for clusters a to d, respectively. The inter-group distances were found to be significantly higher than intra-group SNP distances for each cluster (Fig 5C). Collectively, our data highlighted a common genetic background between commercial strains from ssp aizawai and kurstaki and isolates related to 47 FBOs, corresponding to 18.8% of Bcassociated FBOs (Fig 6). However, insecticides from subspecies israelensis and morrisoni were located out of the four defined clusters a to d, along with two FBO-Bt isolates.

Discussion
Because of a lack of investigations, the involvement of Bt in FBOs is most likely underestimated. To explore this hypothesis, we screened for the presence of Bt within the largest collection of FBO-associated Bc isolates ever analyzed to our knowledge. We report the detection of at least one type of Bt isolate in the context of 49 FBOs, out of 250 Bc-associated FBOs analyzed (19.6%). To address the question of whether Bt could be responsible for the toxic events, the presence of other food pathogens was verified from collected data. In about 39% of Bt-associated FBOs, no other pathogenic microorganism could be detected, except for Bt, making it the most likely causal agent. However, we cannot completely exclude the presence of another microorganism that may not have been found through in vitro detection, but this assumption would concern a relatively high number of toxic episodes (19 FBOs from that point), which seems unlikely. The 49 Bt-associated FBOs affected more than 673 people in total in France between 2007 and 2017, and caused symptoms such as diarrhea and vomiting. Vomiting is commonly considered the result of the Bc emetic toxin, cereulide [47]. Nevertheless, none of the identified Bt isolates were positive for the detection of ces, suggesting the reported vomiting-type symptoms could be correlated with another virulence marker, potentially not described so far. Additionally, almost all Bt isolates harbored genes encoding three major enterotoxins (HblC/D/A, NheA/B/C, and CytK2), predicted to generate diarrhea-type symptoms. Interestingly, Guinebretiere and co-authors showed that the co-occurrence of cytK, nhe and hbl was significantly more often associated with food-poisoning Bc strains than food-related Bc strains [34]. According to previous studies [34], the levels of Nhe and Hbl production of Bt isolates measured in this study are intermediate, when compared with the whole Bc group. However, in another study carried out by the same authors, Bc isolates from phylogenetic group IV-2 (nhe +, cytk+ and hbl+) were found to be associated with elevated cytotoxicity levels on Caco2 cells [48]. This suggests that these enterotoxins are not fully responsible for the cytotoxic damage linked to Bc isolates, and calls again to screen for other virulence markers, previously described or not, along with the cytotoxicity assessment of Bt isolates characterized in this study.
Bt were enumerated from foodstuffs with a median value close to 1x10 3 CFU/g of food, and for about 90% of FBOs, the enumerations were found to be below the alert threshold established by the French Ministry (10 5 CFU/g), concerning bacteria from the whole Bc group. The infective dose of Bc is considered to mainly range from 10 5 to 10 8 CFU/g [49]. Nevertheless, this dose is regularly questioned since it can vary considerably from an outbreak to another one. For instance, low doses of Bc (2.10 2 to 1.10 3 CFU/g) have also been described in foods causing diseases [21,29]. One possible explanation would be that the aggregation of Bc spores strongly interfere with the enumeration of bacilli [21]. Moreover, the infective dose of spores in diarrheal diseases is predicted to be lower than the one of vegetative cells, due to their resistance to the acidic gastric environment [50]. Therefore, our results are consistent with previous findings and suggest that the alert threshold might need to be re-evaluated.
In more than half of cases, the foodstuffs from which Bt were isolated were composed of raw food and in particular tomatoes, whereas non-Bt isolates were 7 to 9 times less frequently associated with raw food. These results are consistent with the Bt prevalence estimations made by Kim et al. who detected Bt in 77% of tested organic vegetables [51]. Similarly, previous studies have demonstrated the presence of high counts (above 10 4 CFU/g) of Bt on fresh vegetables originating from Denmark [52], and some of them were undistinguishable by low resolution genotyping methods from commercial strains such as ABTS-351 [25], suggesting they could be residues of Bt insecticides. Tomatoes and lettuce do not need to be peeled before consumption, contrary to other raw vegetables, such as cucumber. This might explain why these type of food were more frequently incriminated in Bt-FBO, compared to vegetables that are subjected to the same Bt treatments. Taken together, these data reinforce the hypothesis that Bt food contamination may be of agricultural origin. However, we found that in some cases, Bt was reported to be isolated from foodstuffs not subjected to Bt treatments, e.g. Bt isolate 08CEB074, isolated from smoked salmon in FBO 3. The lack of information related to the exact composition of mixed dishes, as well as putative cross contamination between foodstuffs during food processing could easily explain this observation.
In fact, Bt ssp. kurstaki and aizawai are the most widely used bacterial insecticides [53], in particular for the treatment of various crops (including solanaceaous vegetables and cucurbits) and tree fruits [6]. To explore this putative path of contamination, we compared Bt isolated from FBOs to a representative panel of Bt pesticides. All FBO-Bt except for two, appeared to share common features with commercial Bt aizawai and kurstaki strains, including the presence of genes encoding the enterotoxins cytK2, nheA/B/C and hblC/D/A. By contrast, insecticidal Bt ssp. israelensis (BMP144 and AM65-52) and morrisoni (NB-176) are different due to specific genetic signatures, suggesting that they are distinct from FBO-Bt. These results are consistent with the M13-PCR characterization, highlighting similarities of pattern between pesticides belonging to ssp. kurstaki and aizawai, and about 97% of Bt isolates coming from FBOs. Nonetheless, this RAPD-type PCR is a rapid approach for routine comparison of Bc isolates and did not allow us to differentiate commercial Bt strains beyond the subspecies level, leading us to use a genomic approach for further discrimination.
A preliminary phylogeny based on shared k-mers suggested that most FBO-Bt strains we studied were embedded in a major clade of genetically close isolates (phylogenetic group IV), along with a few publicly available genomes of Bt isolated from insects. The isolates from this clade were genetically more distant from isolates from phylogenetic groups VII, then I, VI, V, II and III, successively, which is consistent with previous findings published by A. Bazinet [8].
However, the inclusion of public sequences in the dataset illustrates that the Bt species are diverse and are distributed among four different panC groups (i.e. II, III and V, in addition to IV). This suggests that the panel of FBO-Bt strains analyzed in this study does not reflect the full diversity of environmental Bt isolates, but appears rather to be associated with genetically close subspecies.
Based on an SNP variant analysis, we highlighted that Bt associated with 47 FBOs were specifically organized into four clusters, along with seven pesticide Bt strains belonging to ssp. kurstaki and aizawai. Since the variant calling implemented is restricted to the core genome, we confirmed this distribution with a complementary tree built from accessory genomes. Whereas new tools in WGS analysis have greatly facilitated the investigation on sources of foodborne illnesses, the elaboration of precise guidelines for interpretations remains a complex issue [54,55]. In particular, applying core genome SNP thresholds appears insufficient since the bacterial genome plasticity can result in mutations depending of culture and isolation conditions [56]. Even if authors estimated that pairwise SNP distances between related genomes range from 0 to 4 for Salmonella enterica [54] and Shiga toxin-producing E. coli [57], this statement is not true for long term outbreaks where the outbreak duration and evolution rate induce increasing of pairwise SNP distances between related genomes [55]. Nevertheless, the low computed pairwise SNP distances for Bt isolates within clusters defined in the present study (0 to 2 median pairwise SNP distances) suggested that the isolates were genetically similar and likely originated from a same source. It would be useful for further analysis to investigate the evolutionary capacity of commercial Bt strains, with complementary genomic and laboratory approaches, in order to elaborate interpretation guidelines.
Three pesticide genomes were found to be distinct from FBO-Bt. The corresponding ssp. israelensis and morrisoni are respectively dedicated to control of mosquito larvae and to the protection of potato crops by spraying of the foliage, a part of the plant that is not consumed. Again, the fact that these insecticide strains are distinguishable from FBO-Bt strains is consistent with the hypothesis of an agricultural origin for Bt-associated FBOs. Bc is a ubiquitous and telluric organism; its presence in soils as a natural contaminant is also a possible hypothesis. Our results indicate that Bt was often isolated in the presence of tomatoes, for instance, which are usually cultivated relatively far from the soil. In the hypothetical case of contamination from soil due to raindrop impacts, it seems unlikely that the natural presence of Bt in soils would allow detection of high amounts of bacteria on tomatoes (about 10 3 CFU/g).

Conclusions
This study reported the involvement of Bt in 8 to 20% of Bc-associated FBOs considered for the study, with Bc as the second most common cause of bacterial FBOs in France. Moreover, multiple analyses converged towards the putative agricultural origin of the majority of the identified Bt isolates. Hence, given the high degree of similarity observed between FBO-associated and some commercial Bt isolates, and in the absence of contradictory evidence, we cannot rule out that pesticide Bt may have pathogenic potential. This emphasizes the need for a better understanding of the pathogenicity potential of this agent for non-target organisms, along with the development of specific monitoring tools for Bt traceability in food.  Table. List of phenotypic and genotypic characteristics of 59 representative FBO-Bt isolates and 19 pesticide Bt strains. +/-= detection/no detection of corresponding activity or gene. The phylogenetic groups were assigned according to the partial sequencing of panC [9]. The attribution to M13 groups (named 1 to 3) was established in this study based on similarity Dice coefficients calculated with Bionumerics. (PDF) S3 Table. Nhe absorbance and Hbl dilution index values of a panel of FBO-and commercial Bt. � The Nhe absorbance corresponds to the coloration intensity, obtained with the BDE VIATM kit (3M-Tecra), according to the manufacturer's recommendations. In the tested conditions, the limit of detection corresponds to Nhe Absorbance (414nm) = 0.2. �� The Hbl dilution values correspond to the highest dilution for which Hbl remained detectable, using the immunoenzymatic kit BCET-RPLA (Oxoïd), and according to the manufacturer's recommendations. nd = not detected. Dilution = 1 is the limit of detection of Hbl and corresponds to a detection of 2ng/ml of Hbl components. (PDF) S4 Table. List of 234 isolates used for the k-mer phylogeny. a the attribution to Bacillus thuringiensis (Bt) species was determined by detection of parasporal crystals by phase-contrast microscopy. b the phylogenetic clustering was based on panC sequences similarities [9]. c All the sequencing data obtained from this study are associated with the BioProject PRJNA547495. d,e The strains used for genome assemblies and the SNP calling analysis are indicated with a "x" in the columns "reference assemblies" and "iVARCall2", respectively. (PDF)  [44] and iQtree [43]. The visualization was done using iTOL [46]. (B) Diagrams of the distributions of intra-and inter-clusters log 10 (number of different genes). (PDF) S1 Raw images. (TIF)