Profiling of Bacillus cereus on Canadian grain

Microorganisms that cause foodborne illnesses challenge the food industry; however, environmental studies of these microorganisms on raw grain, prior to food processing, are uncommon. Bacillus cereus sensu lato is a diverse group of bacteria that is common in our everyday environment and occupy a wide array of niches. While some of these bacteria are beneficial to agriculture due to their entomopathogenic properties, others can cause foodborne illness; therefore, characterization of these bacteria is important from both agricultural and food safety standpoints. We performed a survey of wheat and flax grain samples in 2018 (n = 508) and 2017 (n = 636) and discovered that B. cereus was present in the majority of grain samples, as 56.3% and 85.2%, in two years respectively. Whole genome sequencing and comparative genomics of 109 presumptive B. cereus isolates indicates that most of the isolates were closely related and formed two genetically distinct groups. Comparisons to the available genomes of reference strains suggested that the members of these two groups are not closely related to strains previously reported to cause foodborne illness. From the same data set, another, genetically more diverse group of B. cereus was inferred, which had varying levels of similarity to previously reported strains that caused disease. Genomic analysis and PCR amplification of genes linked to toxin production indicated that most of the isolates carry the genes nheA and hbID, while other toxin genes and gene clusters, such as ces, were infrequent. This report of B. cereus on grain from Canada is the first of its kind and demonstrates the value of surveillance of bacteria naturally associated with raw agricultural commodities such as cereal grain and oilseeds.


Introduction
Bacillus cereus sensu lato (B. cereus) is a group of Gram-positive, rod-shaped, spore forming bacteria, which occupy diverse lifestyles and ecological niches, which is facilitate by its diversity at the gene and genome level [1]. These bacteria and their propagules are present in soil, dust, plants, water, feces, and animal guts [2][3][4]. In agricultural, strains of these bacteria can associate with the rhizosphere to promote plant growth and are used as biocontrols to manage insect pests and microbial diseases in wheat fields [1]. Other strains of the bacteria are becoming an area of concern in the food industry because they are able to cause foodborne illness and their particularly well-suited for contaminating dry grains. Bacillus cereus spores and vegetative forms can be inhabitants of many plants, including cereals and their derivatives [1,6]. The objective of the present study was to perform genomic characterization of B. cereus from wheat and flax grains that are produced annually for human and animal consumption.

Enrichment and isolation of B. cereus from raw grains
In 2017 and 2018, wheat and flax grain samples were received from across Canada. Samples were either obtained from the Harvest Sample Program (HSP) of the Canadian Grain commission or were from cargos that had gone through the grain handling system (i.e. processed at grain elevators and transported by truck and/or rail) and are being prepared for export. Wheat samples from HSP belonged to two wheat classes that are defined according to the Canadian grain registration system, Canadian Western Red Spring (CWRS) and Canadian Eastern Red Spring (CERS) and originated from farms across Canada. Canadian Western Red Spring (CWRS) wheat is the most widely grown wheat in Canada and is produced mostly in the prairies (Manitoba, Saskatchewan, and Alberta), while Canadian Eastern Red Spring (CERS) is from Eastern Canada (Atlantic Canada, Quebec, and Ontario). Cargo samples include CWRS grain samples and accompanying dockage samples (any material intermixed with grain, other than kernels of grain), which we analyzed seperately. Upon receipt and prior to any other handling or processing, grain samples were sub-sampled into sample bags, which were then stored at -30˚C. Grain sub-samples were processed for bacterial culturing and enrichment using a modified protocol revised from Health Canada's Compendium of Analytical Methods for the Microbiological Analysis of Foods (MFLP-52, Nov 2014), which is optimized for the enrichment of diverse bacteria from raw commodities. Briefly, to enrich for bacteria in the samples, 50 g of grain was transferred into sterile blender bags with filters (Innovation Diagnostics) and 125 mL buffered peptone water was added to each sample. Samples were then incubated for 24 h at 37˚C. Next, 125 mL of 10% of tryptone soya broth (TSB) was added to the sample bags and samples were homogenized using a pulsifier (Microgen Bioproduct). After the samples were pulsified, an additional 200 mL of 10% TSB was added to each sample, and samples were incubated at 42˚C for 4 h. After incubation, 5 mL of vancomycin (10 μg/mL) and cefsulodin (3 μg/mL) were added to each of the samples, which were then incubated at 42˚C for 16 h. Sample bags were then sealed and thoroughly mixed. The liquid portion of the enriched samples were then aliquoted for bacterial testing.
Enriched samples from 2017 were plated onto Petri-dishes containing Chromagar (Alere Inc, CA) and incubated at 42˚C for approximately 16 h. On Chromagar, B. cereus colonies are expected to appear blue with a white halo, allowing them to be differentiated from other bacteria that may be present. Colonies were then preserved in 15% glycerol at -140˚C. At a later date, 10μl of the frozen culture was re-plated on new plate of Chromagar and confirmed the colony morphology and purity again prior to DNA extraction.
Samples from 2018 were tested for B. cereus directly by real-time PCR. DNA was extracted from the enriched cultures using the Blood and Tissue kit (Qiagen), which was automated using a QIAcube (Qiagen). DNA concentrations were determined using the Qubit dsDNA (double-stranded DNA) Broad Range (BR) assay kit (Invitrogen) and DNA quality was measured spectrophotmetrically by Nanodrop 2000 (Thermofisher Scientific). Samples were tested for the presence of B. cereus by real-time PCR using the Ba primer set (S1 Table) as in [33], with slight modifications. Each 20 μL reaction contained 10 μL SsoAdvanced Universal SYBR green supermix (BioRad), 2 μL of the template DNA, and target specific primers (S1 Table).
Amplification was carried out on a CFX96 Optics module (BioRad) and data were analyzed using CFX Manager software (Bio-Rad).

DNA isolation and PCR testing for toxin genes
Further DNA testing and genome sequencing was performed on 109 strains of B. cereus isolated from the grain samples (S2 Table). Ninety-nine (99) strains of B. cereus originated from Canada Western Red Spring (CWRS) wheat (Triticum aestivum), two (2) strains were from Canada Eastern Red Spring (CERS) wheat, three (3) of them were from wheat dockage, and three (3) were isolated from flax (Linum usitatissimum). Further, two (2) of the previously lab confirmed B. cereus strains (00-0028 and 05-0322) from the National Microbiology Laboratory NML, Winnipeg were also included in the analysis.
Prior to DNA isolation, B. cereus strains were grown and sub-cultured on B. cereus Chromagar plates as described above. To prepare colonies for DNA isolation, 3 mL 0.9% NaCl solution was applied onto the surface of each agar plate and the resulting suspension was aspirated with a pipette and transferred into a reaction tube. Genomic DNA (gDNA) was purified from the isolates using Qiagen Blood and Tissue kit (Qiagen), which was automated using the QIAcube (Qiagen), according to the protocol recommended by the manufacturer for the Gram positive bacteria. All bacterial samples (109) were screened for the presence of five virulence genes (hbID, nheA, ces, cytk1 and cytk2) using real-time PCR as in [34][35][36] (S1 Table), with slight modifications, as described for the Ba primer set.

Whole genome sequencing, assembly, and annotation
To capture the genomic diversity of the suspected B. cereus isolates, we performed whole genome sequencing (S2 Table). Genomic DNA was used to construct NEBNext Ultra II libraries for 94 CWRS samples obtained from harvest year 2018 according to the manufacturer's instructions. These samples were sequenced using a single lane of an Illumina HiSeqX instrument, which generated 150 bp paired-end reads. Paired-end sequencing reads were then trimmed using the Trimmomatic [38] with options to remove Illumina adapter sequences and to have an average quality score of 30. De novo genome assemblies were then generated for each strain using shovill (https://github.com/tseemann/shovill). Two (2) of the suspected cereulide positive samples (BC35N and BC88N) from 2018 were particularly interesting due to their toxin profiles and were suspected to contain the ces toxin genes. These two samples and were instead sequenced using Oxford Nanopore Technology R9 flow-cells, and their genomes assembled using canu [39]. Genomes of the rest of the 13 samples from harvest year 2017 and from the National Microbiology Laboratory (NML), Canada, were sequenced on a PacBio RS II sequencer and their genomes assembled at Genome Quebec, Canada. Assembly statistics were generated for each of the genomes using QUAST [40], and gene predictions and annotations were performed using PROKKA [41]. Multilocus sequence typing (MLST) was then performed using PubMLST database and default paramaters of mlst v2.19.0.

Whole genome comparisons
The genomes were compared to bacterial genome databases using Refseq_masher matches (https://github.com/phac-nml/refseq_masher). The top five matches by Refseq_masher were then downloaded from RefSeq database at the National Center for Biotechnology Information (NCBI) and aligned to each of the genomes using 'nucmer' from the software MUMmer [42], with minimum match -l set to 250. The percent genome alignment to each of the references were then extracted using 'dnadiff' and representative reference genomes were used for clustering analysis in R using the 'ComplexHeatmap' package. Nucleotide variants were identified

PLOS ONE
Profiling of B. cereus on Canadian grain from the 'nucmer' alignments using the genome of the emetic B. cereus strain AH187 (Accession: GCF_000021225.1) as a reference. Variants were filtered to remove variants with a minor allele frequency <0.05. BC14 and BC89 were also removed from the analysis because they were off-target species that had poor alignment and variant calling to the other genomes. Variants were then subject to principal component and 1,000 iterations of hierarchical clustering analysis in R using 'prcomp' and 'pvclust', respectively. Based on the scree plot of the nucleotide variant data, seven components were selected for clustering analysis (S1 Fig).

Analysis of toxin genes, plasmids and antimicrobial genes
To examine the presence/absence of gene homologues with implications in food-borne illness and insecticidal activity (S3 Table), sequences were queried against the whole genome sequences from this study using tBLASTn. BLAST hits were filtered to have >70% alignment length, >50% amino acid sequence identity and <10 −5 e-value and binary matrix was used to assign 1 for the presence and 0 for the absence of the gene. Since cytK1 and cytK2 are similar in sequence, BLASTn was performed using the primer sequences that distinguish the two toxins (S1 Table). In parallel, we mapped reads directly to the genes of interest using Kma v1.3.23 (https://bitbucket.org/genomicepidemiology/kma/src/master/) using the default parameters recommended for each sequencing technology; genes were considered present if identified with a percent identity and coverage of 90% or greater. To detect plasmids and antimicrobial resistance genes, the staramr tool (https://github.com/phac-nml/staramr) was used to scan bacterial genome contigs against the ResFinder, PointFinder, and PlasmidFinder databases.

B. cereus is prevalent on wheat grain
The incidence of B. cereus in grain was assessed by screening flax grains, wheat grains, and dockage from wheat grain samples. The Chromagar culture-based method was implemented in 2017, whereas real-time PCR was used in 2018 (Table 1).
Of the 636 samples tested in 2017, 85.22% tested positive for B. cereus using the Chromagar method ( Table 1). The flax samples from 2017 had the highest level of suspected B. cereus samples; however, the positive rate in Canadian Western Red Spring (CWRS) and Canadian Eastern Red Spring (CERS) samples also had high incidence. Real-time PCR of 508 bacterial samples from 2018 identified fewer positive samples, 56.29% than in 2017 (Table 1). In 2018, the highest incidence of B.cereus positive samples were detected in CWRS HSP and CERS HSP respectively. Nevertheless, we observed the same trend of reduced B. cereus positive samples in the cargo samples, when compared to HSP samples in both years (Table 1). It is unclear if the differences observed between 2017 and 2018 are the due to sensitivity and accuracy differences between real-time PCR and Chromagar detection methods or the result of variation in bacterial levels in the grain from the two years. Nevertheless, our findings indicate that incidence of B. cereus on wheat and flax grains is high, which is expected given the ubiquity of these bacteria in our environment [2]. Our results are consistent with reports from other studies that detected B. cereus in 81% of incoming wheat samples to be used for milling [11]. Is important to note that the results from this study only reflect incidence (presence/ absence) and therefore do not reflect the amount of B. cereus in the sample, since culturing enriched the bacteria present. Future studies of raw samples without culturing/enrichment would be required to provide a true count of the bacteria in the original grain samples. Other studies indicate that, while present, the levels of B. cereus in grain are low, and the number of samples testing positive is reduced after grain conditioning [11]. Our observation of reduced incidence of B. cereus in cargo samples is interesting and suggests that levels of the bacteria may be decreasing as grain continues through the supply chain. While B. cereus has been reported in cereals previously [6,11], it is possible that bacterial levels may be higher on the grain immediately after harvest due to sources of bacteria from the field, and that these levels decrease over time as a result of safe storage, handling, transport, and cleaning of the grain.

Grain harbors diverse strains of B. cereus
To better capture the diversity of the suspected B. cereus isolates from this study, we performed whole genome sequencing of 109 samples, which were sub-cultured on Chromagar medium (S2 Table). These samples were selected from both 2017 and 2018 harvest years. Genomes were then de novo assembled, which generated a median genome size of~6.2Mb and GC content of~35% (S4 Table).
To identify candidate species for each sample, genomes were compared to available genomes from the RefSeq database at NCBI using RefSeq_masher and from MLST analysis (S5 and S6 Tables). The majority of the samples were identified to be similar to species within the B. cereus group by RefSeq_masher and were identified to be from the B. cereus scheme by MLST (S5 and S6 Tables). To validate the species assignments and better identify the relationship between the samples and those from publically available genomes, we downloaded the genomes for the top RefSeq_masher hits and performed whole genome alignments to each of our newly assembled genomes by MUMmer (Fig 1; S7 Table). Whole genome alignments confirmed that that the majority of the samples had high genome alignment to those from B. cereus group sequenced from other studies. Most of the samples fell within two closely related groups, containing 40 (Group 1) and 43 (Group 2) samples, respectively. The remainder of the samples (Group 3) were much more diverse and were more distantly related to Groups 1 and 2 (Fig 1). The sample relationships were in agreement with a principal components analysis of nucleotide variants, which demonstrated a close relationship between Groups 1 and 2, with some separation by Principal Component 3 (Fig 2). Group 3 was much more diverse and contained several subgroups (Figs 1 and 2).
Three of the samples (BC35N, BC88N, and BC108P) had over 97% genome alignment to the reference genomes for NC7401 and AH187, which are emetic strains of B. cereus and were isolated from cases of foodborne illness [43]. Four strains (BC93, BC52, BC30, and BC77) had ~70% genome alignment to B. wiedmannii strain BAG6X1, and formed a subgroup within Group 3. For some of the other samples, the species assignment was less clear. For example, BC69 aligned well to Group 1 and 2; however, alignment was too low for it to be clearly placed in these species, with the greatest alignment (~65%), being to strains G9842 and Lr7/2 (Fig 2C  blue arrow; S7 Table). Curiously, the genome size was quite large at 8.9Mb, but the GC content was 38%. Likewise, six strains (BC87, BC43, BC01, B76, BC 79, and BC67) all had reduced alignment to available references; however, unlike BC69, these strains had poor alignment to most Group 1 and 2 references, but had some alignment to Group 3 references (Figs 1 and 2).
A small set of the samples were identified by RefSeq_masher to be bacteria from other genera or species, including Enterococcus, Actinobacter, Brevibaccilus, and other Bacillus species (S5-S7 Tables). The presence of other bacterial species in the samples was further supported by differences in genome size and GC content when compare to those that were identified to be from the B. cereus group, suggesting that multiple bacteria may be present in some samples (S4 Table). Curiously, two of the samples appear to be pure samples from more distant bacteria, including one isolate related to Brevibacillus sp (BC89) and one isolate of B. megaterium (BC14), which formed an outgroup in our genome alignment analyses (Fig 1) and were  Table).

PLOS ONE
Profiling of B. cereus on Canadian grain excluded from the variant analysis (Fig 2; S5-S7 Tables). Both of these organisms have previously been reported in the soil/rhizosphere and contribute to complex microbe-microbe and/ or microbe-plant interactions, with Brevibacillus sp. having biocontrol and bioremediation potential and B. megaterium promoting plant growth [44,45]. While detection and sequencing of these off-target organisms was both interesting and unexpected, the majority of the samples that were sequenced contained diverse strains of B. cereus (Figs 1 and 2).

Analysis of genome and gene content
B. cereus is genetically diverse and known for its presence/absence variation in genes involved in disease, toxin production, and antimicrobial resistance (AMR) [33][34][35][36]. This variation is partly due to the presence of plasmids, which are very common in B. cereus; some B. thuringiensis strains may carry more than 15 plasmids [46]. In addition, plasmid size varies greatly and the plasmid profile of the different strains does not always match their phylogeny. The B. cereus genome assemblies were analysed for plasmids using PlasmidFinder for replicon typing/subtyping (S9 Table). PlasmidFinder supported the presence of plasmids in the genomes, including pBMB67 in BC26, which was previously identified in B. thuringiensis and was implicated in cell-cell signaling and regulation of cellular processes [47].
Some B. cereus strains are known to produce enterotoxin and emetic toxins that have implications in disease; as such, we performed real-time PCR, BLAST, and Kma analyses of some toxin genes. B. cereus toxins that are associated with foodborne illness include hemolysin BL (partially encoded by hbLD), the non-hemolytic enterotoxin Nhe (partially encoded by nheA), cytotoxin K (cytk), and cereulide (ces) (S10 Table) [34][35][36]. The agreement in the PCR and tBLASTn results were 90.8% for hblD, 97.2% for nheA, 84.4% for cytk1, 67.0% for cytk2, and 100% for ces. This was similar to the analysis comparing the PCR and Kma results, which were in 89.0, 98.2, 87.2 and 67.9 and 100% agreement with the PCR results, respectively. The high concordance in the results for hblD, nheA, and ces demonstrate the reliability of PCR-based methods in the detection of most of these genes (S10 Table). The reduced agreement in the cytk results may be because of additional diversity in the cytk alleles. Our real-time PCR results showed that a large proportion of the B. cereus group strains from grain carry the nheA (96.33%) and hbID (94.5%) genes, while 12.8%, 42.2%, and 2.8% strains carried cytk1, cytk2, and ces respectively (Fig 3). The frequent occurrence of nheA and of hblD genes was comparable with the previous studies [48,49]. The cytK gene was less frequent than nhe, and hbl, which is also in agreement with previous studies of food (37%) and diarrheal type of food poisoning (73%) [48,49]. Curiously, cytK2 was less common in samples from Group 1, despite their genomic similarity to Group 2 (Fig 3). The ces gene, encoding emetic toxin and known as cereulide synthetase, was detected in only three samples of CWRS, one from 2017 and two from 2018 (BC35N, BC88N, and BC108P). These are the three sample samples that clustered together within Group 3 and had high genome alignment to the emetic strains NC7401 and AH187 (Figs 1-3). This finding indicates that emetic toxin producing isolates of B. cereus on grain are rare.
To investigate coincidence of toxins, all strains were divided into 11 different toxin genes profiles according to the presence or absence of toxin genes ( Table 2). Pattern I (51 isolates, 46.36%) was the most common pattern and contains enterotoxin genes (nheA and hblD), which was followed by pattern II (35 isolates, 32.11%), which carried both nheA and hblD plus cytk2 genes. This finding is comparable to other studies that identified that the majority of strains (79 of 147; 53.7%) encoded the nheAB and hblDA genes. Patterns V, VI, VIII, X and XI were rare and each pattern consisted only one isolate (0.91%). Patterns II, III, IV, V, X and XI were positive for 3 or more enterotoxin genes, including nheA, hblD, ces or cytKs. Interestingly, the cereulide positive isolates belonged to three different toxin patterns (Patterns V, VI and X), despite their close relationships.
Our analysis of the genes involved in the production of insecticidal crystals (Cry genes; S3 Table) indicate that none of the samples contained the genes investigated. These genes are important for the biopesticide activity of B. thuringiensis in commercially available products registered for use in Canada [50]. The absence of these genes in the samples from this study suggest that the isolates we recovered are not those from commercial products used in agriculture to control insect pests.
In addition to toxin genes, AMR genes are also important because they could impact the fitness of the bacteria and the treatment strategy in cases of disease. The ResFinder tool was used to identify several AMR genes belonging to different classes of antimicrobials. The majority of

PLOS ONE
Profiling of B. cereus on Canadian grain B. cereus samples (61.5%, 67/109) carried either fosB1 or fosB1gene, which confers resistance to fosfomycin (Fig 3; S10 Table). Interestingly, all of the Group 2 samples carried fosfomycin resistance genes, while only a subset of Group 1 samples, and none of the Group 3 samples carried the resistance genes (Fig 3). Genes responsible for resistance to vancomycin (VanC1XY, VanC3XY, VanC3XY and VanC2XY) were detected in fifteen isolates, which were broadly distributed across the different B. cereus groups (Fig 3). Further, five isolates were positive for ISA (A), which is responsible for lincomycin, clindamycin, dalfopristin, pristinamycin 11, virginiamycin, and quinupristin. Tetracycline resistance encoded by tet(S), and the aac(6')-Iid gene that confers resistance to gentamycin, were present in only three isolates each. Taken together, while fosfomycin resistance was common, particularly in samples from Group 2.
The plasmids, toxin genes, and AMR genes identified in this study provide important insights into the status of B. cereus on grain. However, it is important to note that vancomycin and cefsulodin were used as part of our bacterial enrichment method. Therefore, samples may not proportionately represent the B. cereus that may be occurring on grain. There may also be additional genes beyond those reported here that provide resistance to these antimicrobial compounds. In addition, contamination by other bacteria, including Enterococcus and Acinetobacter species, may have resulted in misidentified features in some samples (S5-S7 Tables). Nevertheless, our findings indicate that some toxin genes (nheA and hblD) and AMR (fosfomycin) are common in the B. cereus samples on grain, while others, such as ces are less frequent.

Conclusions
Our study is amongst only a few studies that have investigated B. cereus on raw grain. While we observed that these bacteria were abundant on grain in 2017 and 2018, we noticed a reduction in incidence in later stages of the supply chain. Isolation and whole genome sequencing revealed that these bacteria are diverse and belong to several different species. Further genomic characterization revealed that some of these isolates have the potential to produce toxins; however, strains that produce cereulide are rare. Some of these bacteria may also have resistance to fosfomycin. Since B. cereus is abundant in our natural environments and can occupy different niches, including roles in promoting plant growth and health, it will be important to continue to describe and monitor this group of bacteria and unravel its diversity.