Devon W Kavanaugh

Benjamin Glasset

Rozenn Dervyn

Cyprien Guérin

Sandra Plancade

Sabine Herbin

Anne Brisabois

Pierre Nicolas

Nalini Ramarao
email: nalini.ramarao@inrae.fr

Nalini Rama Rao

New genetic biomarkers to differentiate non-pathogenic from clinically relevant Bacillus cereus strains

Keywords: Bacillus cereus, pathogenicity, genetic biomarkers

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Bacillus cereus is the third causative agent of food-borne-outbreaks (FBO) in Europe [START_REF] Glasset | Large-scale survey of Bacillus cereus-induced food-borne outbreaks: epidemiologic and genetic characterization EuroSurveillance[END_REF]. B. cereus can induce two types of gastrointestinal diseases, leading to generally mild and self-limiting emetic or diarrhoeal syndromes, although several cases of severe infections have been reported [START_REF] Fagerlund | Toxin production in a rare and genetically remote cluster of strains of the Bacillus cereus group[END_REF]. B. cereus also induces systemic infections leading to patient death in approximately 10% of cases [START_REF] Bottone | Bacillus cereus, a volatile human pathogen[END_REF][START_REF] Ramarao | Two unrelated episodes of Bacillus cereus bacteremia in a neonatal intensive care unit[END_REF][START_REF] Gaur | Bacillus cereus bacteremia and meningitis in immunocompromised children[END_REF][START_REF] Lotte | Virulence Analysis of Bacillus cereus Isolated after Death of Preterm Neonates, Nice, France, 2013[END_REF][START_REF] Chan | Infective endophthalmitis caused by Bacillus cereus after cataract extraction surgery[END_REF]. B. cereus is also a source of central nervous system infections and other systemic infections especially in newborns [START_REF] Bottone | Bacillus cereus, a volatile human pathogen[END_REF][START_REF] Cormontagne | Bacillus cereus Induces Severe Infections in Preterm Neonates: Implication at the Hospital and Human Milk Bank Level[END_REF]. Recent epidemiological studies show that the number of cases of serious B. cereus infections is largely underestimated [START_REF] Glasset | Bacillus cereus, a serious cause of nosocomial infections: epidemiologic and genetic survey[END_REF]. The pathogenic potential of B. cereus is extremely variable, with some strains being harmless and others lethal.

B. cereus possesses several toxin genes, such as nhe, hbl a n d c y t K [START_REF] Fagerlund | Toxin production in a rare and genetically remote cluster of strains of the Bacillus cereus group[END_REF][START_REF] Ramarao | The pore-forming haemolysins of Bacillus cereus: a review[END_REF]. These toxins provide an indication of the strain toxicity potential but are not sufficient, alone, to discriminate hazardous from harmless strains [START_REF] Glasset | Bacillus cereus, a serious cause of nosocomial infections: epidemiologic and genetic survey[END_REF][START_REF] Guinebretière | Enterotoxigenic profiles of food-poisoning and food-borne Bacillus cereus strains[END_REF][START_REF] Martinez-Blanch | Development of a real-time PCR assay for detection and quantification of enterotoxigenic members of Bacillus cereus group in food samples[END_REF][START_REF] Ramarao | Advanced Methods for Detection of Bacillus cereus and Its Pathogenic Factors[END_REF]. Indeed, several studies have shown that Nhe production by hazardous strains is variable and that non-pathogenic strains can also produce it in large quantities [START_REF] Glasset | Large-scale survey of Bacillus cereus-induced food-borne outbreaks: epidemiologic and genetic characterization EuroSurveillance[END_REF][START_REF] Martinez-Blanch | Development of a real-time PCR assay for detection and quantification of enterotoxigenic members of Bacillus cereus group in food samples[END_REF]. Moreover, these toxins do not appear to be suitable markers for strains causing non-gastrointestinal infections [START_REF] Glasset | Bacillus cereus, a serious cause of nosocomial infections: epidemiologic and genetic survey[END_REF]. B. cereus produces other toxins such as haemolysin II (HlyII), the metalloproteases InhA1, InhA2 and the cell wall peptidase FM (CwpFM), which may also be involved in pathogenicity [START_REF] Tran | Structural Modeling of Cell Wall Peptidase CwpFM (EntFM) Reveals Distinct Intrinsically Disordered Extensions Specific to Pathogenic Bacillus cereus Strains[END_REF][START_REF] Tran | Bacillus cereus immune escape: a journey within macrophages[END_REF][START_REF] Tran | Hemolysin II is a Bacillus cereus virulence factor that induces apoptosis of macrophages[END_REF][START_REF] Cadot | InhA1, NprA and HlyII as candidates to differentiate pathogenic from non-pathogenic Bacillus cereus strains[END_REF][START_REF] Haydar | InhA1-Mediated Cleavage of the Metalloprotease NprA Allows Bacillus cereus to Escape From Macrophages[END_REF]. The emetic form of B. cereus food poisoning is caused by the peptide cereulide [START_REF] Ehling-Schulz | Identification of emetic toxin producing Bacillus cereus strains by a novel molecular assay[END_REF], which represent less than 1% of the FBO strains of B. cereus [START_REF] Glasset | Large-scale survey of Bacillus cereus-induced food-borne outbreaks: epidemiologic and genetic characterization EuroSurveillance[END_REF][START_REF] Ehling-Schulz | Identification of emetic toxin producing Bacillus cereus strains by a novel molecular assay[END_REF][START_REF] Hoton | The cereulide genetic determinants of emetic Bacillus cereus are plasmid-borne[END_REF].

To date, the above described determinants were not sufficient to completely explain the virulence of B. cereus [START_REF] Glasset | The cytotoxic potential of Bacillus cereus strains of various origins[END_REF] and there are currently no markers, either alone or in combination, to differentiate pathogenic from non-pathogenic strains. In this work, we took advantage of a well characterized collection of 95 B. cereus strains and compared pathogenic (FBO and clinical) with non-pathogenic strains. We identified a combination of four as yet undescribed biomarkers, wherein their presence/absence allows an accurate identification of clinical B. cereus strains.

These findings constitute a huge step in the understanding of the B. cereus pathogenic potential and complexity and may provide tools to better assess the risks associated with B. cereus contamination.

Methods

Isolate information

This study includes 39 B. cereus strains associated with foodborne illness [START_REF] Glasset | Large-scale survey of Bacillus cereus-induced food-borne outbreaks: epidemiologic and genetic characterization EuroSurveillance[END_REF], 35 strains isolated from human patients following systemic or local infections [START_REF] Glasset | Bacillus cereus, a serious cause of nosocomial infections: epidemiologic and genetic survey[END_REF] and 21 non-pathogenic strains [START_REF] Guinebretière | Enterotoxigenic profiles of food-poisoning and food-borne Bacillus cereus strains[END_REF][START_REF] Kamar | Pathogenic Potential of Bacillus cereus Strains as Revealed by Phenotypic Analysis[END_REF] (Sup Table 1). We have previously shown a correlation between cytotoxicity and virulence [START_REF] Glasset | The cytotoxic potential of Bacillus cereus strains of various origins[END_REF]. Nevertheless, although these strains had previously been shown to be weakly cytotoxic to human cells and to have reduced virulence in an insect infection model, this does not rule out their potential ability to produce symptoms in specific vulnerable populations.

RNA extraction

The transcriptome study by RNAseq was carried out on 15 strains representative of the three collections (Sup Table 2) in triplicates. Bacterial cultures were incubated in BHI medium at 30°C in microaerophilic condition (5% O2-15% CO2-80% N2) at pH 7 until entry into stationary growth phase. Samples were centrifuged at 12,000 g for 3 min at 4°C and placed immediately at -80°C until processing. The bacterial pellets were re-suspended with 200 μl of 10 mM Tris-HCl at pH 8 + 4 μl of lysozyme at 50 mg/ml and incubated at 37°C. Total RNA was extracted with the HPRNA kit (High Pure RNA Isolation Kit; Roche) as previously described [START_REF] Porrini | Implication of a Key Region of Six Bacillus cereus Genes Involved in Siroheme Synthesis, Nitrite Reductase Production and Iron Cluster Repair in the Bacterial Response to Nitric Oxide Stress International[END_REF]. The RNA integrity was measured by the RIN (RNA Integrity Number) and were between 7 and 10. The mRNA were enriched with the RiboZero Kit (Illumina). The sequencing of the mRNA was carried out by the I2BC platform (CNRS, Gif-sur-Yvette). Directional and paired libraries were prepared with the Illumina scriptseq kit and the sequencing was performed on an Illumina Nextseq machine.

Transcriptome sequencing analysis

Sequencing quality was assessed using FastQC, and adapter sequences and low-quality base pairs were removed using cutadapt (version 1.9) [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. Reads were further trimmed in 3′ using sickle (version 1.33, option "-x" and default values for all other parameters, implying a Phred quality cutoff of 20). In absence of whole genome sequences for the 15 strains, the cleaned reads were mapped against a repertoire of allelic variants for 23,815 genes aiming at accounting for the pangenome of B. cereus group. This repertoire was obtained by single-linkage clustering based on the results of an all-against-all blastn comparison (version 2.2.26, e-value cut-off 1e-5) [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] of 519,931 CDSs extracted from the 91 annotated complete genomes available at the time of analysis for B. cereus group in Genbank. Pairs of CDSs that aligned over at least 70% of the length of the shortest sequence and with at least 75% nucleotide sequence identity were grouped in the same cluster, which resulted in 23,815 clusters representing distinct genes. Reads were mapped using bowtie2 (version 2.2.6, options "-N 1 -L 16 -R 4") [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] whose results were converted to bam format using SAMtools version 1.9 [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF]. Read counts on each allelic variant were obtained using HTSeq-count (version 0.6.1) [START_REF] Anders | HTSeq--a Python framework to work with high-throughput sequencing data[END_REF] and summed over allelic variants to obtain a single read count per gene per sample. To cope with sequence similarity between allelic variants of a same gene and fragmentation of the reference according to gene boundaries, R1 and R2 reads were aligned independently and use of HTSeq-count option "-a 0" allowed to count reads that aligned equally well on several allelic variants of a same gene. Of note, since bowtie2 mapped each read on a single allelic variant, reads could not be counted more than once in the sum. Expression levels expressed as log2 scaled rpkm (reads per kilobase per million mapped reads) were produced by the R package "edgeR" (version 3.11) using the mean length of the genes in the cluster and a prior count of 1.

Raw transcriptomic data and differential expression analysis are accessible through GEO Series accession number GSE168681

(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171128).

Statistical model

The strategy for statistical analysis of RNAseq data was to select genes to predict whether a strain is pathogenic y=1 or not y=0 and evaluate the prediction accuracy. We considered the logistic regression model with lasso penalty implemented in the R-package "glmnet", which allows the selection of a limited subset of genes whose expression is associated with strain pathogenicity [START_REF] Engebretsen | Statistical predictions with glmnet[END_REF]. The package glmnet provides an interval cross validation procedure to select the penalty constant, which determines the number of selected genes.

The prediction accuracy of the procedure was evaluated in a cross-validation framework where splitting in training and validation sets preserves the matching of the three replicates of each strain. For each replicate, the model provides a probability ẑi to be pathogenic, and we considered the average value over the three replicates as the prediction probability of the strain. The predicted pathogenicity status is set to zero if the prediction probability is smaller than 0.5 and 1 otherwise.

Biomarker screen by PCR

The 7 marker genes were retrieved from at least 20 sequenced B. cereus strains from NCBI databases and aligned by CLC Main workbench7 software to identify two regions conserved across the strains. Within these regions, 20 bp primers were designed using the Beacon Designer software. For the majority of the selected genes there were no perfectly conserved sequence and some bases had to be replaced with R (A/T), Y (C/T) or W (A/T) for primer design (Sup Table 3).

For all the strains of the collection, a single colony was picked, resuspended in 100 µL Tris-EDTA NaCl buffer (TEN) and incubated at 98°C for 10 min. After centrifugation, 1 µl of supernatant was used as DNA matrix. The PCR mixture contained 1 µl DNA matrix, 0.5 µM primer (forward and reverse), 10 µL DreamTaq Green PCR Master Mix (2X) (Thermo Scientific) in a final volume of 20 µL. PCR fragment sizes were revealed on 1.5% agarose gels containing Midori Green, and visualised by a UV imaging device.

AUC analysis to select combinations of biomarkers

The PCR data were pooled into a presence (1) /absence (0) table, which was then used as input for ROC-AUC analysis facilitated by the web-based suite of tools hosted at www.combiroc.eu. The ROC-AUC analysis determines the combination of biomarkers, which will best differentiate the classes of samples input ('disease' versus 'non-disease' groups). Sets of biomarkers were selected based on their performance in sensitivity or specificity alone, or in combination as the AUC metric.

Potential hits were filtered at 85% specificity and 85% sensitivity.

Results

RNAseq analysis

We obtained between 9-15 million reads per samples with 90% correctly paired. The overall alignment rate was over 85%. The analysis enabled the creation of a read counts table based on gene expression levels for each sample (Figure 1). The dispersion of the sample count values was homogeneous and the biological triplicates clustered well together. We identified 3276 genes in the core transcriptome, which represents approximately 65% of the genes in each strain.

Identification of 7 biomarkers by logistic regression analysis

A Mann-Whitney-Wilcoxon nonparametric rank test with a classical 5% of qvalue did not allow the prediction of significant differences in gene expression among the strain collections (not shown). Thus, to identify markers that could potentially differentiate pathogenic from nonpathogenic strains, we performed a penalized conditional logistic regression with the lasso method on the entire counting table to select relevant genes for the prediction of pathogenic potential. By applying the prediction model to the 11,179 genes with the selected penalty constant of 0.01, only 7 genes were selected (Table 1).

With the RPKM values of these 7 genes (Sup Table 4), a prediction in a cross-validation framework among the 15 strains, leads to 13 well classified strains (estimated probability ẑi value below 0.5 for non-pathogenic and above 0.5 for pathogenic strains) and two misclassified strains, one false positive (NP strain PF predicted as pathogenic) and one false negative (pathogenic FBO strain 12CEB01BAC predicted as NP) (Table 2).

Validation of the biomarkers on a large strain collection

Initially, for the first 15 strains, the presence of the 7 selected genes was further assed by PCR (Table 3). These data revealed that when a gene showed no expression by transcriptomic analysis, the gene was actually absent from the strain. Thus, the identification of these 7 biomarkers was based on gene presence/absence, rather than mRNA expression. As such, an approach centred on gene detection was chosen for the screening of the large bacterial collection with the 7 genes selected (Table 3) and to determine the area under the curve (AUC), specificity, and sensitivity of possible combinations of the selected biomarkers.

1-FBO vs NP

For the FBO strains, the best combination of biomarkers able to differentiate NP from FBO strains was obtained with 4 biomarkers (Figure 2A). With this combination, the best AUC was 0.768, the sensitivity 0.69 and the specificity 0.773. Therefore, we obtained some false positive (NP strains that appear pathogenic), and some false negative (FBO strains that appear NP). Taken together, the general trend for the FBO identification was an overall low AUC among the tested combinations, thus preventing their accurate differentiation.

Nevertheless, we identified that several FBO strains were lacking almost all biomarkers. These FBO strains primarily belong to the phylogeny group IV (table 3). We thus performed an additional AUC analysis after the removal of all strains of the phylogeny group IV of the collection (FBO and NP). The results were significantly improved and the best combination resulted in an AUC above 0.9 and with significantly improved sensitivity or improved specificity. But a combination resulting in sensitivity and specificity above 0.9 was not determined (Figure 2B).

2-NP vs clinical strains

Regarding the clinical strains, the best results were achieved with a combination of 4 biomarkers with an AUC of 0.955, sensitivity of 0.9 and specificity of 0.86. Therefore, the analysis concludes that an accurate differentiation between clinical and non-pathogenic strains can be obtained by using these biomarkers (Figure 2C). These two combinations allowed the accurate discrimination between the two strain populations. Some markers have the same occurrence within the strain collection [START_REF] Gaur | Bacillus cereus bacteremia and meningitis in immunocompromised children[END_REF][START_REF] Lotte | Virulence Analysis of Bacillus cereus Isolated after Death of Preterm Neonates, Nice, France, 2013[END_REF][START_REF] Chan | Infective endophthalmitis caused by Bacillus cereus after cataract extraction surgery[END_REF] and were therefore interchangeable during the AUC analysis. Thus, the best combinations of biomarkers are: 1, 2, 3, 5 (or 6 or 7). The genes are named, adhB, agrC, thiJ, BCQ_PI180 (or gshAB or BCQ_PI181).

As a conclusion, a suitable combination of 4 biomarkers has been found to create a robust and accurate test to differentiate clinical from non-pathogenic strains, with an AUC of 0.955, given that test results above 0.9 are considered excellent.

Discussion

The emergence of B. cereus as a foodborne pathogen and as an opportunistic pathogen has intensified the need to distinguish strains of public health concern. The pathogenic potential of B. cereus is extremely variable, with some strains being harmless and others lethal. Currently, due to the lack of validated and standardized analytical methods, only the presence of B. cereus is usually investigated in foods or clinical samples at a species-level. Over the years, new methods have been developed with the leading principle to detect and distinguish B. cereus from others Bacillus group members by a time-saving and in-situ analysis [START_REF] Manzano | Molecular methods to evaluate biodiversity in Bacillus cereus and Bacillus thuringiensis strains from different origins[END_REF], genotyping using highresolution melting analysis [START_REF] Antolinos | Development of a highresolution melting-based approach for efficient differentiation among Bacillus cereus group isolates[END_REF], the use of multi-locus sequence (MLST) [START_REF] Didelot | Evolution of pathogenicity in the Bacillus cereus group[END_REF] or the classification of the strains according to their affiliation to a phylogenetic group that offers a first useful indicator of risk [START_REF] Guinebretière | Enterotoxigenic profiles of food-poisoning and food-borne Bacillus cereus strains[END_REF]. Nevertheless, MLST analysis of the 53 strain sequences included in this study revealed that 21% belonged to the sequence type ST26, and approximately 11% to an undetermined ST (not shown), while >40% of the strains were identified as belonging to PanC clade III (Table 3). As such, the ST types and PanC classifications were unable to completely explain the grouping of the strains.

Here, we report new markers characteristic of pathogenic B. cereus strains, which detection requires only PCR, and is thus independently of growth conditions. We could indeed show that the simple presence/absence of the gene was as discriminant as its expression value by transcriptomic analysis. We further calculated the AUC, specificity and sensitivity obtained using the combination of these 4 biomarkers to discriminate between our large B. cereus collection inducing various pathologies. CombiROC results demonstrate that clinical strains were more efficiently separated from the non-pathogenic strains than the FBO strains.

Regarding the FBO strains, to improve the analysis, strains belonging to the phylogenetic group IV were removed, thus allowing a significant improvement in strain differentiation. This might prove very useful for food industries to better communicate the risks of B. cereus food contamination and to take the appropriate measures for decontamination while preventing or minimizing economic loss. Nevertheless, this implies a two step-test with a first panC phylogenetic attribution followed by a biomarker test.

By contrast, regarding the clinical strains, the combination of 4 biomarkers allowed the identification of a strong differentiation test with an AUC of 0.955, sensitivity of 0.9, and specificity of 0.86. Thus, a global test with a strong AUC (above 0.9) and increased sensitivity (rare false negative) could be proposed to accurately discriminate between clinical and harmless strains. As such, our new findings may be relevant to gain additional knowledge on the strains found in hospitals and healthcare settings. Table 1. List of 7 selected biomarkers with gene position (on the reference genome pAH187_270 -NC_011655.1) and putative function.

Legends of figures and tables

Table 2. Estimated probability ẑi for the 15 strains. A logistic regression model with lasso penalty was applied to select the penalty constant, which determines the number of selected genes.

Then prediction accuracy of the procedure was evaluated in a cross-validation framework. For each replicate, the model provides a probability ẑi to be pathogenic, and we considered the average value over the three replicates as the prediction probability of the strain. The predicted non-pathogenicity corresponds to a ẑi smaller than 0.5 and the predicted pathogenicity corresponds to ẑi above 0.5.

Table 3.

Presence/absence of biomarkers among non-pathogenic (green), FBO (blue) and clinical (beige) strains. The presence of each biomarker gene was assessed by PCR in all strain of the collection. If the gene was present, a score of 1 was attributed (green boxes), if the gene is absent, a score of 0 is attributed (red boxes). Table 2: Estimated probability ẑi for the 15 strains. A logistic regression model with lasso penalty was applied to select the penalty constant, which determines the number of selected genes.

Then prediction accuracy of the procedure was evaluated in a cross-validation framework. For each replicate, the model provides a probability ẑi to be pathogenic, and we considered the average value over the three replicates as the prediction probability of the strain. The predicted non-pathogenicity corresponds to a ẑi smaller than 0.5 and the predicted pathogenicity corresponds to ẑi above 0.5.

NP

Figure 1 .

 1 Figure 1. RNAseq heatmap. Heatmap representation of expression levels (log2 rpkm) across the pangenomic repertoire of 23,815 genes (rows) and the 45 samples (columns). Dendrograms are built by hierarchical clustering with average-link. The 3,272 genes with signal in all strains are indicated by grey bars. Non-pathogenic strains are indicated in black and pathogenic strains in red.

Figure 2 .

 2 Figure 2. CombiROC analysis results. The presence/absence matrix resulting from PCR detection of biomarker sequences was analyzed by CombiROC. (A) Foodborne outbreak strains (FBO) versus non-pathogenic; (B) FBO versus non-pathogenic strains, excluding phylogenetic group IV. Links best sensitivity performance, right highest specificity; (C) clinical versus non-pathogenic strains.

Table 1 :

 1 list of the 7 selected biomarkers with gene position and putative function

	, araC, gshAB	0.768 0.692 0.773

Acknowledgments

We want to thank Sandrine Auger and Mylène Sperry for their help in the primer design. We thank Valentin Loux for bioinformatics support.

COI statement

The authors declare no conflict of interest.

Funding

This work was supported by the European EJP Toxdetect project from the European Union's Horizon 2020 research and innovation program under Grant Agreement No 773830 and by the Comue Paris Saclay Idex Program n°CDE-2018-002323 -IRE 2018-0021.

Author Contributions

DK, BG, RD: performed experiments, analyzed data, manuscript writing; CG, SP, PN: analyzed data; SH, AB: supervision; NR: initial concept, supervision, analyzed data, writing of manuscript, funding sources.

Table 3. Presence/absence of biomarkers among non-pathogenic (green), FBO (blue) and clinical (beige) strains. The presence of each biomarker gene was assessed by PCR in all strain of the collection.

If the gene was present, a score of 1 was attributed (green boxes), if the gene is absent, a score of 0 is attributed (red boxes).