Evaluating the Safety of Bacillus cereus GW-01 Obtained from Sheep Rumen Chyme

Bacillus cereus is responsible for 1.4–12% food poisoning outbreaks worldwide. The safety concerns associated with the applications of B. cereus in health and medicine have been controversial due to its dual role as a pathogen for foodborne diseases and a probiotic in humans and animals. In this study, the pathogenicity of B. cereus GW-01 was assessed by comparative genomic, and transcriptome analysis. Phylogenetic analysis based on a single-copy gene showed clustering of the strain GW-01, and 54 B. cereus strains from the NCBI were classified into six major groups (I–VI), which were then associated with the source region and sequence types (STs). Transcriptome results indicated that the expression of most genes related with toxins secretion in GW-01 was downregulated compared to that in the lag phase. Overall, these findings suggest that GW-01 is not directly associated with pathogenic Bacillus cereus and highlight an insightful strategy for assessing the safety of novel B. cereus strains.


Introduction
Bacillus cereus is a ubiquitous Gram-positive, facultative anaerobic, spore-forming, rod-shaped bacterium [1,2], which has been isolated from various sources, such as soil, water, sewage sludge, air, human skin, food, and vegetation [3].The safety concerns associated with B. cereus have been controversial, as the bacterium acts as both a pathogen for foodborne diseases [4] and probiotic in humans and animals [5][6][7].

RNA Sequencing and Transcriptomics Analysis
GW-01 was inoculated in LB medium at 30 • C and 150 rpm, after which cells at lag, logarithmic, and stationary phases were harvested via centrifugation (12,000× g, 10 min, 4 • C) for RNA extraction.As previously reported [35], the RNA extraction and sequencing transcriptomics were performed by Shanghai Personal Biotechnology Co., Ltd.(Shanghai, China).In this process, RNA was extracted and purified using the Trizol Reagent (Invitrogen, Carlsbad, CA, USA) and the RNeasy mini kit (Qiagen, Venlo, The Netherlands).The purity and integrity were assessed using the NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA) and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively.Bacterial ribosomal RNA (rRNA) was depleted using the ribo-Zero rRNA removal kit (Illumina, San Diego, CA, USA), and RNA library was prepared with Illumina TruSeq Stranded total RNA kit.Further, paired end sequencing was performed on Illumina HiSeq X Ten (China Shanghai personalbio Co., Ltd., Shanghai, China), and transcriptomic data were aligned and mapped onto the complete genome of the GW-01 [35] strain using bowtie2 (https: //bowtie-bio.sourceforge.net/index.shtmlaccessed on 7 July 2024).Differential gene expression analysis was normalized to reads per kilobases per million reads (RPKM) values, and RNA-Seq output data were analyzed for statistical significance using the Rockhopper software (http://cs.wellesley.edu/~btjaden/Rockhopperaccessed on 7 July 2024).

Data Analyses
The data from toxicity assessment assay were analyzed by using the SPSS 16.0 software (SPSS, Chicago, IL, USA), presented as the mean ± standard deviation (SD).One-way ANOVA was carried out to detect the statistical significance.Statistical significance was set at p < 0.05.The sequence alignment was performed employing Mega-X (MEGA 1.0).

Phylogenetic Analysis of Bacillus cereus GW-01
The reputation of B. cereus as a significant foodborne disease pathogen has resulted in the extensive study of the population structure and phylogenetic relationships of the isolates within the group have using diverse typing methods, to help monitor the evolution of strains or to identify clones responsible for disease or food poisoning outbreaks [46][47][48][49][50][51].Thus, it is necessary to study the safety of GW-01 through phylogenetic analysis.Interestingly, phylogenetic analysis showed that it was difficult to distinguish the Bacillus cereus from Bacillus thuringiensis, Bacillus anthracis, Bacillus weihenstephanensis, Bacillus mycoides, Bacillus pseudomycoides, Bacillus gaemokensis, Bacillus manliponensis, Bacillus cytotoxicus, Bacillus toyonensis, Bacillus bingmayongensis, and Bacillus wiedmanni [52], and most research groups have proposed them as Bacillus cereus sensu lato (s.l.).The pathogenicity in this group of bacteria is very diverse, thus making it challenging to distinguish bacterial species that cause disease or food poisoning outbreaks [33,53].In this study, the Average Nucleotide Identity (ANI) value between strain GW-01 and ATCC 14579 was 91.64%, and the ANI value of these 55 strains ranged from 90% to 98% (Supplemental Document S1).Based on computational biology, some of these 55 strains did not belong to same species although all of them were classified under B. cereus in the NCBI by classical microbial taxonomy.Several research groups have proposed the classification by computational biology based on whole-genome sequences of Bacillus cereus sensu lato (s.l.).[46,54].However, most food scholars still identified the strain by classical microbial taxonomy, which was the main reason for choosing Bacillus cereus species with complete sequence information from the NCBI for further analysis.
In this study, the GC content (%) and genome size of the 55 B. cereus strains ranged from 34.89% to 35.55%, and 5.16 to 5.56 Mb, respectively (Figure 1a,b and Table S2).GW-01 has the highest GC content (35.55%), and its genome size is 5.24 Mb, suggesting that its evolution was different from that of other strains.Its genome size makes it difficult to analyze its pathogenicity, due to the notion that well-adapted pathogenic bacteria generally contain smaller genomes (due to less variable selection pressure) than environmental isolates [55,56].However, this notion is debatable.Moreover, bacterial pathogenicity is associated with the habitats and hosts [57], and there is no apparent relationship between the evolution, habitats, and hosts of these 55 strains (Figure 1a,b).
The pathogenicity of B. cereus is reportedly connected to their evolution [61,[66][67][68], and multiple strains from diverse habitats and hosts show genomic variations in terms of gene content, while the essential functions for a specific species are conserved as core genes that are shared among different strains [57].In this study, phylogenetic analysis based on the single-copy core genes and the core genes was performed to determine the relationship of GW-01 and ATCC 14579 with other 53 known strains of B. cereus spp.(Figure 1a,b).These analyses led to clustering of the 55 B. cereus strains into six major groups (I-VI).Zervas et al. [69] have also suggested that Bacillus cereus sensu lato (s.l.) strains could be classified into six major groups by PanC and multiple loci sequence typing (MLST) characterization, and ANI-based phylogenetic trees.Our pan-genomic analyses revealed the presence of 30,868 genes in 55 B. cereus strains; 2221 and 22,874 of these constituted the core and accessory genomes, respectively (Figures S1 and S2).As shown in Figure 1a, most strains in cluster II and V were isolated from USA, while cluster IV and VI strains were mainly found in China and Republic of Korea (East Asia).A phylogenetic analysis based on the core genes was performed to describe the relationship of GW-01 and ATCC 14579 with the other 53 B. cereus strains (Figure 1b).Cluster I-V strains were mainly found in China and Republic of Korea (East Asia).These results indicated geographical association with evolution of B. cereus spp.However, the phylogenetic clusters based on core genes of the 55 B. cereus strains were not identical to those that phylogenetically divided the strains based on single-copy genes.GW-01 clustered together with CC-1 and branched near the base of the CC-1 and M3 (Figure 1b), which is consistent with the phylogenetic analysis based on single-copy gene (Figure 1a).This suggests that GW-01 likely arose before and after CC-1 and M3, respectively.Interestingly, Torres Manno, Repizo, Magni, Dunlap and Espariz [54] suggested that M3 and CC-1 were assigned to Bacillus pacificus and Bacillus paranthracis, respectively.Different species of Bacillus cereus sensu lato (s.l.) locating in same cluster was also reported by other studies [53,61,70].GW-01 was not found to cluster with the known toxic strains of B. cereus spp.except for FORC_021 and FORC_005.Thus, not all the bacterial species within a given genetic group exhibit the same level of pathogenicity, and the distinction between pathogenic and innocuous strains remain unclear, for the entire B. cereus group.The pathogenicity of GW-01, CC-1, and M3 remain uncertain, but they were not clustered together with strain 03BB87, FORC_005, FORC_021, S2-8, ATCC 14579, and FM1, which have been previously reported as pathogenic.

Prophage, CRISPR, T3SS Effector Protein, VFDB, CARD, and Genomic Island Predictions
During evolution, the presence of prophages is an important driving factor responsible for the expression of pathogenic properties.The bacteriophage-mediated DNA transfer can convert a non-pathogenic strain into a pathogenic variety through prophage-encoded toxins, surface alterations, or increasing resistance to human immunity [71][72][73].Marraffini Luciano and Sontheimer Erik [74] have reported the predicted role of CRISPR in limiting phage insertions.Thus, both prophage and CRISPR play an important role in virulence evolution in bacterial pathogens.However, prophage and CRISPR system in these 55 strains were not reported as related with the bacterial pathogens in this study (Supplemental Document S2).The highest GC content of prophage (66.3%) and CRISPR (68.1%) were both appeared in GW-01, indicating their difference with other strains.Moreover, T3SS effector protein, VFDB, and CARD gene were analyzed (Supplemental Document S3), but no significant association was found with the bacterial pathogens.
Genomic islands (GIs) contribute to the evolution of pathogenic microorganisms [75].Therefore, to understand the pathogenic attributes of the 55 B. cereus strains, predictions of genomics islands (Gis) were performed using IslandPath-DIOMB, SIGI-HMM, IslandPicker, and IslandViewer.Highest number of genomic islands (18) were present in B. cereus MLY1 and DLOU-Tangshan, and least number of genomic islands (5) were identified in B. cereus AR156 (Table S3 and Supplemental Document S4).Notably, in B. cereus GW-01, five genomic islands were identified.Moreover, the highest antibiotic resistance numbers were 7 in B. cereus G1-1.The highest and least pathogen-host interaction (PHI) numbers were 15 in B. cereus 03BB87 and DLOU-Changhai each, and 2 in B. cereus ISSFR-9F, respectively.In strain GW-01, gene, victors gene, antibiotic resistance, and PHI numbers were 200, 12, 2, and 8, respectively.For pathogenic bacteria, the PHI number in strain 03BB87, FORC_005, FORC_021, S2-8, and FM1 are 15, 7, 9, 8, and 8, respectively; the genomic island number in these strains were 13, 11, 14, 14, and 9, respectively.Desvaux et al. [76] suggested that the functional features of PHI played an important role in bacteria pathogenicity.The gene description is shown in Supplemental Document S1.Strain GW-01 does not exhibit some of the PHIs known to be possessed by pathogenic bacteria.It suggested that strain GW-01 pathogenicity was different with these pathogenic strains.

Secretion Pathway of Toxin in GW-01 Based on Transcriptome Analysis
Although our results indicated that the evolutionary property of GW-01 was different from pathogenic B. cereus, its pathway of toxin secretion needed further investigation.The transcriptome results of GW-01 at lag, logarithmic, and stationary phases are shown in Figure 2. The differently expressed genes (DGEs) between lag, logarithmic, and/or stationary phases are shown in Figure 2a, Figure 2b, and Figure 2c, respectively.Compared to the lag phase, upregulated and downregulated genes were 1033 and 933 at the logarithmic phase, and were 1753 and 1768 at stationary phase, respectively.Upregulated and downregulated genes were 611 and 1289 at stationary phase compared to the logarithmic phase.The terms related with B. cereus growth and proliferation significantly changed as growth phase prolonged (Figure 2d-f).Interestingly, salmonella infection and flagellar assembly terms between the lag and logarithmic phases were significantly changed, suggesting that GW-01 infection might associated with its growth.Duport et al. [77] indicated that B. cereus exoproteome was mainly produced during the logarithmic phase.
In Bacillus cereus sensu lato (s.l.) strains, the diarrheal symptoms are caused by poreforming toxins hemolysin BL (HBL), non-hemolytic enterotoxin (nhe), cytotoxin K (CytK), and hemolysin I (CLO), as well as phospholipases, proteases, and the entD protein [33].Almost all B. cereus strains contain at least one type of toxic protein [4,26,78], with nhe and HBL being present in approximately 100% and 50% of all enteropathogenic B. cereus strains, respectively [4,[79][80][81].In this study, only nhe (A, B, and C) was found to be present in the strain GW-01.Subsequently, we attempted to analyze bacterial pathogenicity by performing the sequence alignment of toxicity-associated genes.However, it was still difficult to distinguish the pathogenicity, as the toxicity-associated genes in all the strains of B. cereus spp.were similar.Our analysis thus revealed that the similarity of nhe between strain GW-01 and 54 B. cereus strains was higher than 94% (Table S4).Toxins from certain pathogenic bacteria must cross the cell membrane(s) to gain access to their site of action in the target host cell [82].Thus, protein secretion is important in the pathogenicity of B. cereus strains.Toxic protein secretion is accomplished by flagellar export apparatus (FEA), resulting its deficient secretion in non-flagellated strains [83,84].In this study, flagellar assembly terms in logarithmic, and stationary phases were significantly downregulated compared to the lag phase.At least six integral proteins (flhA, flhB, fliO, fliP, fliQ, and fliR in Salmonella) have been reported to participate in the flagellar export [85,86], which were substantially homologous with the type III secretion system (T3SS) proteins [87,88].In this study, T3SS analysis in 54 B. cereus strains revealed that the distributions of T3SS proteins between strain GW-01 and the known pathogenic strains are not regular (Table S5A).Furthermore, the amount of secreted HBL in a B. cereus mutant lacking flhF was shown to be significantly reduced, for which a signal recognition particle (SRP)-like GTPase is involved in the regulation of the number and arrangement of flagella on the bacterial surface [89].The nhe secretion in B. cereus is influenced by T3SS proteins, flhF and flhA [90].Similarity of genetic determinants in B. cereus strain GW-01 with respect to these genes in the rest of the 54 strains was higher than 90% (Table S5A).The expression of flhF, flhA, and other genes related with T3SS at the lag phase was downregulated compared to those at logarithmic, and stationary phases (Table S6), suggesting the toxins secretion in GW-01 through T3SS did not increase as growth phase prolonged.In Bacillus cereus sensu lato (s.l.) strains, the diarrheal symptoms are caused by po forming toxins hemolysin BL (HBL), non-hemolytic enterotoxin (nhe), cytotoxin K (Cy and hemolysin I (CLO), as well as phospholipases, proteases, and the entD protein [ Almost all B. cereus strains contain at least one type of toxic protein [4,26,78], with nhe HBL being present in approximately 100% and 50% of all enteropathogenic B. cer strains, respectively [4,[79][80][81].In this study, only nhe (A, B, and C) was found to be pres in the strain GW-01.Subsequently, we attempted to analyze bacterial pathogenicity It also has been demonstrated that HBL, nhe, and cytK secretion is dependent on the general secretory (Sec) pathway [82], yet genes encoding hbl and cytK were not detected in the strain GW-01.The proteins related with Sec path include secB, secA, srp, yidC, secYEG, ftsY, secDF-yajC, pmf, and spaces [91].Interestingly, secB and secA were not found in GW-01, and the expression of other genes were downregulated (Table S7).These results indicated that the toxins in GW-01 might not have been secreted in Sec pathway.Huang et al. [92] have reviewed that the genetic determinants, spo0A, comER, plcR, codY, abrB, rpoN (Sigma 54), and flhA (Flagella) are associated with toxin production in B. cereus spp.Thus, it is necessary to investigate the similarity of these genetic determinants between strain GW-01 and 54 B. cereus strains in order to understand the pathogenicity of B. cereus strain GW-01.
The probable relationship between B. cereus toxins and regulated genes was hypothesized (Figure 3) to understand the pathogenicity of B. cereus strain GW-01.The spo0A protein regulating either the spore formation or biofilm development pathways in B. cereus [93], is responsible for the contamination of food during manufacture and processing [92].Moreover, biofilm formation and early sporulation in B. cereus were also positively associated with the comER [94].In this study, the similarity of spo0A and comER between the strain GW-01 and the rest of the 54 B. cereus strains was found to be higher than 97% and 94% (Table S5B), respectively, suggesting that spo0A and comER have a highly conserved amino acid sequence.Lindbäck et al. [95] have reported that codY, a global regulator, may indirectly upregulate enterotoxin production in B. cereus, while, flhA an important factor for toxic proteins secretion in B. cereus, encodes a component of flagellum-apparatus formation [92].Our sequence alignment studies have revealed that the similarities of codY and flhA between the strain GW-01 and the rest of the 54 B. cereus strains were higher than 99% and 97%, respectively, indicating their highly conserved amino acid sequence (Table S5A,B).Furthermore, both the distribution and similarity of sinI and sinR in all the 53 B. cereus strains, as well as the strain GW-01, were irregular (Table S5A,B).These results indicated that spo0A, comER, codY and flhA are highly conserved while sinI and sinR failed to analyze the pathogenicity in the B. cereus strain GW-01.Hayrapetyan et al. [96] have reported that rpoN (Sigma 54) downregulates the toxin (nheA) production and sinR expression in B. cereus ATCC 14579.In this study, only strain GW-01 contained rpoN (Sigma 54) genes, indicating that nhe production in B. cereus strain GW-01 might be limited.However, plcR, which activates the expression of hbl, nhe, and cytK [92], was also present in the strain GW-01 (Figure 3).Interestingly, all the 52 B. cereus strains except for the strain NJ-W contain plcR genes, while the similarity between the GW-01 and NJ-W strains was found to be only 12.82%.Further, most known pathogenic B. cereus strains possess the papR genes; their similarities with strain GW-01 were determined to be higher than 91% (Table S5B).However, Hayrapetyan, Tempelaars, Nierop Groot and Abee [96] has reported that the expression of plcR and its regulon members was downregulated in the rpoN mutant.Thus, we speculated that the rpoN (Sigma 54) gene was associated with the toxin production either by regulating by plcR/papR or through unknown mechanisms.In this study, the expressions of genes related with toxin synthesis were downregulated except the expression of SinI and SinR at stationary phase compared with the lag phase, and the logarithmic phase, suggesting that the expression of Nhe was downregulated (Figure 3).Heilkenbrinker et al. [97] have suggested that toxicity would be reduced when the ratios of nheB and nheC are lower and higher than 50:1 and 5:1, respectively.In this study, the ratios of nheB and nheC expression at lag, logarithmic, and stationary phases were 8.7:1, 21.7:1, and 3.4:1, respectively (Table S8).Thus, we concluded that the pathogenicity of GW-01 at the lag and logarithmic phases might be higher than that at stationary phase.

Conclusions
B. cereus is a well-known foodborne pathogen and the toxin produced by it has been the cause of food poisoning outbreaks worldwide.However, certain strains are widely used as

Figure 2 .
Figure 2. Transcriptome analysis of GW-01 at lag, logarithmic, and stationary phases.Volcano (a) and KEGG enrichment (d) between the lag phase and the logarithmic phase; volcano plot (b) KEGG enrichment (e) between the lag phase and the logarithmic phase; volcano plot (c) and KE enrichment (f) between the lag phase and the logarithmic phase.

Figure 2 .
Figure 2. Transcriptome analysis of GW-01 at lag, logarithmic, and stationary phases.Volcano plot (a) and KEGG enrichment (d) between the lag phase and the logarithmic phase; volcano plot (b) and KEGG enrichment (e) between the lag phase and the logarithmic phase; volcano plot (c) and KEGG enrichment (f) between the lag phase and the logarithmic phase.

Figure 3 .
Figure 3. (a-c) represent the possible pathway and gene expression of toxin synthesis in lag, logarithmic, and stationary phases of GW-01.Green indicates downregulation, and red indicates upregulation.