Genetic changes are introduced by repeated exposure of Salmonella spiked in low water activity and high fat matrix to heat

WGS is used to define if isolates are “in” or “out” of an outbreak and/or microbial root cause investigation. No threshold of genetic differences is fixed and the conclusions on similarity between isolates are mainly based on the knowledge generated from previous outbreak investigations and reported mutation rates. Mutation rates in Salmonella when exposed to food processing conditions are lacking. Thus, in this study, the ability of heat and dry stress to cause genetic changes in two Salmonella serotypes frequently isolated from low moisture foods was investigated. S. enterica serovars S. Agona ATCC 51,957 and S. Mbandaka NCTC 7892 (ATCC 51,958) were repeatedly exposed to heat (90 °C for 5 min) in a low water activity and high fat matrix. No increased fitness of the strains was observed after 10 repeated heat treatments. However, genetic changes were introduced and the number of genetic differences increased with every heat treatment cycle. The genetic changes appeared randomly in the genome and were responsible for a population of diverse isolates with 0 to 28 allelic differences (0 to 38 SNPs) between them. This knowledge is key to interpret WGS results for source tracking investigations as part of a root cause analysis in a contamination event as isolates are exposed to stress conditions.

WGS results interpretation depends on pairwise SNP distances between isolates, bootstrap support of the branches in a phylogenetic tree and the tree topology 15 . Among these three parameters, genetic distances estimated as SNP and in case of cg/wgMLST analysis as allelic differences, often forms the initial basis to determine the similarity between isolates. Interpretation of SNP/ allelic values to conclude on the similarity of isolates can be challenging as absolute thresholds do not exist or it is not feasible to prescribe one. The prevailing knowledge on the genetic distances to determine the similarity between isolates is mostly derived from previous outbreak investigations and reported mutation rates. In a cross-sectional study, Wang et al. 16 investigated 6,351 Salmonella isolates originating from 2,196 facilities and concluded that isolates from within a facility will have fewer (< 20) SNP differences (probability = 0.66), though larger differences could also occur. The outbreak investigation of S. Agona infections among infants in 2017 indicated that the outbreak isolates clustered within a maximum distance of 26 SNPs and these isolates originated from cases with gastrointestinal symptoms between April and December of 2017 17 . Salmonella mutation rates of 9.3 × 10 -8 per nucleotide/year for the accumulation of core SNPS (or 0.44 SNPs per genome/year) calculated for S. Agona 18 and 2.2 × 10 -7 per site/year (or 1.01 SNPs per genome/year) calculated for S. Enteritidis 19 are reported. An extrapolation of these mutation rates does not correlate with the number of SNP differences observed between the outbreak S. Agona isolates within a time period of 8 months, indicating the potential role of other factors. However, there is limited information on the Salmonella evolution rate in the farm to fork continuum which is essential to interpret WGS results of isolates of food origin and its associated environment. Though WGS analysis of longitudinal set of Salmonella isolates recovered from a single facility might provide insights into the mutation rates, such studies are generally lacking and even more it is difficult to link the impact of specific environmental or food process conditions on the genetic changes. In this context, the objective of this study was to evaluate the ability of simulated sublethal food processing conditions to induce genetic changes in Salmonella as to the best of our knowledge this phenomenon has not been studied. Two model serotypes frequently associated with low moisture foods were used in this study. S. Agona has been associated with foods such as aniseed-containing herbal tea 20 , infant milk products 17 and cereals 21 while S. Mbandaka has been associated with cereals and is a commonly found Salmonella serovar 22,23 . The strains were repeatedly exposed to heat in a low moisture matrix to assess its impact to introduce genetic changes. In addition, phenotypic parameters namely, log reduction, lag time and sublethal injury levels were studied to evaluate if changes in fitness could be observed and if so to potentially correlate with the introduced nucleotide changes.

Materials and methods
Experimental design. The overall experimental procedure is outlined in Fig. 1. Previously, we have shown that a maximum of only 1 SNP was introduced after 10 sub-culturing steps in TSA/ Columbia agar (37 °C, 24 ± 2 h) for Salmonella enterica serovars (Tennessee, Hadar, Typhimurium and Enteritidis) 24 . While mutations can be caused by sub-culturing steps, based on our previous findings, this is less likely to happen and thus we premised that mutations observed following this protocol would largely be induced by exposure to stress conditions.  Inoculation of low water activity and high fat matrix (hereinafter matrix). The inoculum was diluted 100-times in TS and 1 mL was inoculated in 200 µL droplets to 100 g of sterile (irradiated Synergy Health, Däniken, CH) dry matrix in a stomacher bag (1% inoculum level to obtain a target level approximately 10 6 CFU/g). Afterwards, the bag was massaged by hand until no clumps were observed anymore. Each time 3 replicates were analyzed after inoculation to ensure the inoculation was done homogeneously. The standard deviations for each of the 10 inoculated matrices was on average 0.14 and 0.17 log cfu/g for S. Agona and S. Mbandaka, respectively, showing the inoculation was homogenous. The matrix consisted of dried digested animal byproducts, an ingredient used for the coating of dry pet food, with a composition of 3.1% moisture, 60% protein and 25% fat content. The a w of the non-inocuated matrix was 0.208. The inoculated matrix was hermetically closed in an aluminum bag (Vacopack, Jehud, Israel) and stabilized in an incubator at 25 °C during 24 h following which the a w increased to 0.275.
Heat treatment of the spiked matrix. Thermal cells as described by Rachon et al. 25 were used to carry out the heat treatments. For each replicate, 1.5 g of inoculated matrix was transferred to a thermal cell (0.8 mm deep,48 mm diameter). Three replicates were prepared. The control thermal cell, containing 1.5 g non-inoculated matrix, had an incorporated built-in platinum thermocouples (Pt 100) designed and supplied by Nestlé Research (Lausanne, Switzerland). In each trial, temperature profiling was conducted and the core temperature of samples were recorded using a data logger (PicoLog TC-08; St Neots, UK). The 3 replicates with inoculated matrix and the control thermal cells were immersed in a water bath at 90 °C for 5 min. After the heat treatment, the thermal cells were returned to room temperature using a water bath (25 °C, 1 min).

Enumeration on selective (XLD) and non-selective (TSA) medium. From the thermal cells, one
g of heat-treated matrix was transferred to a stomacher bag and tenfold serial dilutions were prepared in TS to carry out the enumeration on TSA and Xylose Lysine Deoxycholate (XLD) (ThermoScientific Oxoid, Hampshire, UK) agar (24 h incubation at 37 °C). The difference in counts on TSA before and subsequent to the heat treatment were used to calculate the log reduction caused by each heat treatment cycle. The counts (CFU/g) on XLD, and TSA were used to calculate sublethal injury levels after each heat treatment cycle using the following formula: [(Count on TSA) -(Count on XLD)]/ (Count on TSA). Enumeration of the inoculum level (S. Agona and S. Mbandaka) was also performed on TSA and XLD before it was used to spike the matrix.

Determination of lag time and inoculum preparation for 10 consecutive heat treatment cycles.
A growth curve of the inoculum was determined by adding 5 mL of the 10 -7 dilution from the inoculum (approximately 10 10 -10 11 CFU/mL) to 45 mL of preheated (37 °C) BHI broth. Enumeration was carried out on TSA at 0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8 and 24 h to determine the growth curves and to calculate the lag phase of the inoculum using ComBase (A Web Resource for Quantitative and Predictive Food Microbiology. University of Tasmania; USDA Agricultural Research Service. https:// data. nal. usda. gov/ datas et/ comba se-web-resou rcequant itati ve-and-predi ctive-food-micro biolo gy).
Following the heat treatment, 5 mL of tenfold dilution from each sample was added to 45 mL of preheated (37 °C) BHI. Enumeration was carried out on TSA at 0, 0.5, 1, 2, 3, 4, 5, 6, 7, 8 and 24 h to determine growth curves and to estimate the lag phase by ComBase. The inoculum for the next round of heat treatment was prepared by growing a lawn culture of the 10 -2 and 10 -3 dilutions from a single inoculated sample on TSA (37 °C, 24 h) and subsequently following the procedure described above. This new inoculum was used to spike a new batch of matrix for the next cycle of heat treatment cycle. In total, 10 consecutive heat treatment cycles were performed on the two strains of Salmonellae. For each heat treatment cycle, 3 colonies picked from TSA and XLD, where possible, were sequenced from the inoculum and the 3 replicates of the heat treated matrix.

Statistical analysis.
A global p-value (ANOVA) as well as pairwise comparisons (post-hoc comparison without comparisons for multiplicity of tests) was calculated with R version 3.6.1 (R: A language and environment for statistical computing. R Core Team. R Foundation for Statistical Computing, Vienna, Austria,2019, ISBN 3-900,051-07-0, URL http:// www.R-proje ct. org/) for the log reduction on TSA and the lag time.
WGS analysis. Three individual colonies from TSA and XLD respectively isolated from the inoculum and heat-treated samples were inoculated into 4 mL of BHI broth and incubated at 37 °C for 6 h. After incubation, an aliquot (one mL) of BHI was centrifuged at 5000 × g, 5 min. The pellet was store at -20 °C until the DNA extraction was performed according to Gimonet et al. 26 and sequenced with HiSeq as described by Portmann et al. 24  For S. Mbandaka, at least 3 colonies were present on XLD for all cycles but for S. Agona for cycles 5, 6, 7 and 10, less colonies were taken since not always 3 colonies for each replicate were available and no colonies were taken for cycles 1 to 4.
BioNumerics v7.6.3 (https:// www. appli ed-maths. com/ bionu merics) was used to carry out the wgMLST analysis. Assembly-free and assembly-based allele calling were used to calculate the allelic differences using the default settings. The Salmonella scheme consists of 15,874 loci. A minimum spanning tree (MST) was created to illustrate the relatedness between all sequenced isolates for S. Agona and S. Mbandaka.
High quality SNP (hqSNP) pipeline developed by the Center for Food Safety and Applied Nutrition (CFSAN) at the U.S. Food and Drug Administration (CFSAN SNP Pipeline v.1.0.1) was used for SNP calling on the isolates 27 .
Maximum-likelihood phylogenetic tree were built with GARLI (Version 2.01.1067 28 ) on the SNP analysis results.
Identification of genes in which genetic changes were observed for S. Agona and S. Mbandaka. Basic Local Alignment Search Tool (Blastn version 2.4.0 + 29 ) was used to identify the genomic regions containing a SNP comparing the concerned nucleotide region with the ATCC strains. Accession numbers of the ATCC strains are NZ_CP019183 for S. Mbandaka and NZ_AOZX01000005, NZ_AOZX01000007 NZ_ AOZX01000011, NZ_AOZX010000636, NZ_AOZX01000039, NZ_AOZX01000043 and NZ_AOZX01000064 for S. Agona. SNPs were investigated that occurred at least 10 times along all S. Agona and S. Mbandaka isolates respectively.
Visualization of the genomic regions was done using Artemis Comparison Tool (ACT) 30 .

Results
Phenotypic characteristics following exposure of S. Agona and S. Mbandaka to heat and dry stress. The reduction in counts (CFU/g) by repeated heat treatments (90 °C for 5 min) in an inoculated low water activity, high fat matrix are presented in Fig. 2. Between 0.2 and 2.1 log reduction was obtained in the case of S. Agona. The 3 rd heat treatment cycle resulted in a markedly lower log reduction compared to the other cycles. In contrast, the 10 th cycle caused a significantly higher log reduction compared to the 9 previous heat treatment cycles. Overall, the log reduction of S. Mbandaka was similar (P > 0.05) among the 10 heat treatment cycles, except for cycle 4 where no log reduction was observed. The reduction levels were consistent for S. Mbandaka after the 4 th heat treatment cycle on TSA but 1.5 ± 0.4 log reduction was observed on the selective medium XLD (data not shown). The average % of sublethal injury levels before and after heat treatment are mentioned in Table 1. Sublethal injury levels varied in the inoculated matrix (before heat treatment) for both Salmonella strains. After heat treatment, > 98% of the surviving cells were sublethally injured and the variability between replicates was significantly lower.
Additionally, the lag time was calculated for (i) the inoculum that was used to spike the matrix and (ii) subsequent to the heat treatment of the spiked matrix. The lag time of the inoculum was stable along the 10 cycles (Figs. 3 and 4) and was significantly lower (P < 0.05) in comparison with the lag times of both strains after the heat treatment. The lag time of S. Agona after heat treatment was significantly lower (P < 0.05) after cycle 8 and 9 compared to the other cycles (Fig. 3). The lag time of S. Mbandaka from the first cycle was significantly higher (P < 0.05) compared to the 9 subsequent cycles (Fig. 4).  Figs. 5 and 6, respectively. The isolates associated with the same treatment cycle are marked with the same color in both Figs. 5 and 6. The isolates from the initial cycle form the core of the MST and ML tree. As evidenced from the figures, the differences grew with increasing heat treatment cycles. The isolates did not group    www.nature.com/scientificreports/ per cycle e.g. the isolates from the 10 th cycle are dispersed throughout the MST. A maximum of 28 allelic or 38 SNP differences were observed between 160 S. Agona isolates that were sequenced and analyzed. In Table 2, the maximum allelic and SNP differences noticed per cycle are mentioned. The number of differences accrued with every additional heat treatment cycle performed thus showing an increasing trend per cycle. The MST obtained from wgMLST analysis of S. Mbandaka is shown in Figs. 7 and 8 represents the ML tree generated from SNP analysis. The maximum number of allelic and SNP differences is 8 and 19, respectively, among the 210 sequences obtained from S. Mbandaka isolates compared to 28 for S. Agona. The number of identical isolates (zero allelic/ SNP differences), represented by the core circle in Figs. 7 and 8, is higher when compared to S. Agona. Similar to S. Agona, an increasing trend of allelic differences was observed with the  www.nature.com/scientificreports/ increasing number of treatment cycles, though it was less pronounced. A maximum of 19 SNP differences was noticed at cycle 7, and it decreased to 10 differences by cycle 10. Using wgMLST analysis, a maximum of 8 allelic differences was obtained after cycle 7 and it did not further increase. Majority of the SNPs were unique in nature (91% for both serovars, 201/220 for S. Mbandaka and 722/794 SNPs for S. Agona) i.e. they were observed only once, suggesting that most of the mutations were most likely random in nature. In order to gain a preliminary understanding about the biological relevance of these SNPs, we investigated the SNPs that occurred minimum 10 times following exposure to stress conditions. Thus, for S. Mbandaka 5 SNP positions and for S. Agona 8 were selected for further investigation. The base change position and impact on the amino acid coded along with the gene or intergenic region implicated is presented in the supplementary Table S1. Four non-synonymous changes and one SNP in the intergenic region were observed for S. Mbandaka, whereas, S. Agona had 2 SNPs in the intergenic regions, 1 synonymous and 5 non-synonymous amino acid changes. A total of 2 out of 4 and 2 out of 5 SNPs creating non-synonymous changes in S. Mbandaka and S. Agona, respectively, were noticed in genes coding/ controlling flagellar functions.

Discussion
Subtyping tools provide evidence or can bring new clues to direct the epidemiological investigation in an outbreak situation or for the root cause analysis of a pathogen contamination event in a food plant. These analytical results help to determine which isolates are "in" or "out" of the scope of an investigation. Traditional subtyping tools such as PFGE are being increasingly replaced by WGS as it provides a higher resolution to discriminate between pathogen isolates. In 2017, a S. Agona outbreak caused by contaminated infant milk products among infants was identified in France 17 . The confirmed cases were analyzed by WGS and the SNP analysis revealed that the outbreak isolates clustered within a maximum of 26 SNPs 17 . In a different study from Germany, the potential similarity between a local S. Agona isolate from a feed sample to the French outbreak isolates was investigated 31 . Results indicated that the feed sample isolates differed by at least 40 SNPs in comparison to the French outbreak isolates and thus they were considered unrelated. However, the results also showed that an outbreak cluster linked to tea (2003) and a different cluster of coconut (1994) differed by only 21 SNPs to the isolate from the feed sample. It has been shown that fixed genetic cut-off values cannot be assigned due to the difference in evolution forces between organisms, the environmental exposure history of the organism and the context (e.g. duration of the outbreak, ability of an environment to support growth impacting generation time) 12,15 . In addition, the organisms can also be exposed to conditions causing sublethal stress injury during food processing that may induce genetic changes likely due to selection pressure. As low moisture and heat are two stress conditions that are present in a low moisture food producing facility, we questioned whether the exposure of Salmonella, known to be associated with low moisture foods, to these conditions could induce genetic changes. In this study, we therefore investigated the effect of heat stress on two Salmonella strains S. Agona and S. Mbandaka inoculated in a high fat, low water activity matrix.
The effect of the subsequent heat treatments upon the fitness of the strains was evaluated by the phenotypic characteristics such as inactivation by log reduction, lag time and sublethal injury. We observed a maximum of 2.1 log reduction in the matrix consisting of 25% fat and a water activity of less than 0.3. It has been shown that Salmonella is rapidly killed in liquid conditions but in matrices with high fat content and low water activity, the inactivation by heat is reduced 32 . Each heat treatment cycle resulted in > 98% levels of sublethal injury. The lag time did not show significant changes over the ten repeated heat cycles except for S. Mbandaka where a significant decrease in lag time was observed after one heat treatment. In conclusion, the phenotypic characteristics remained mostly unchanged over 10 repeated heat treatments in a low moisture, high fat matrix. This observation is in contrast to the study of Knöppel et al. 33 , who showed by serial passaging, and thus repeated exposure to, different growth media resulted in increased fitness. After 500-1000 generations in parallel cultures in 4 different media, 83 genetic changes consisting of amino acid substitution, deletion, duplication, frameshift mutation, intergenic mutation, frameshift reversion or pseudo reversion, non-sense mutation and synonymous mutation were identified for S. enterica LT2 by the authors.
On the genotypic level, genetic changes were observed for S. Agona and S. Mbandaka after the exposure to heat treatment in the low water activity, high fat matrix. It has been shown that S. Tennessee, Typhimurium and Enteritidis showed a maximum of 1 SNP difference after 10 sub-culturing events 24 . Maximum 3 substitutions after four passages in S. Montevideo was observed 34 . After 100 subcultures of S. Typhimurium and S. Newport, between 0 and 3 allelic differences were observed by wgMLST 35 . Thus, sub-culturing of different Salmonella serovars induced a maximum of 3 genetic changes. Therefore, we can conclude that significantly more genetic changes were introduced by exposure to repeated heat treatment, as observed in this study, than genetic changes that would have been expected by sub-culturing without stress. Additionally, an increasing number of allelic/ SNP differences was observed after repeated exposure to heat. The number of genetic changes was higher for S. www.nature.com/scientificreports/ Agona compared to S. Mbandaka. The genetic diversity within Salmonella depends upon serovar e.g. S. Newport 36 exhibit a large genetic diversity versus S. Enteritidis which has a low genetic diversity 19 . It should be noted that the laboratory simulations mimic stress conditions that might occur in an environmental production plant but does not reflect the amount of generations expected in these circumstances. After each heat treatment cycle, strains were incubated in non-selective medium for 24 h at 37 °C while in real production areas the chances of such occurrence are unlikely. From the S. Agona MST, it seems that the randomly introduced genetic changes remains in the population and subsequent cycles introduce additional genetic changes in the isolates. As a result, the diversity of sequences increases with an increase in the number of heat treatment cycles. Similar observations, but to a lesser extent, were noticed for S. Mbandaka, i.e. the population with different genomic sequences was higher after cycle 10 compared to cycle 1. Overall, similar patterns were observed with the ML tree obtained for both serovars following the SNP analysis. In terms of genetic differences, SNP differences were almost always slightly higher than allelic distances but were comparable to each other for both serovars in different treatment cycles, except in cycle 1 for S. Agona. The higher number of differences noticed with the SNP analysis can be explained by the fact that about 95% of the genome is normally considered in the analysis, whereas the only coding genes are considered in wgMLST 11 . However, both these methods have been shown to provide concordant results for other organisms [12][13][14] and the findings of our study confirms this observation for Salmonella enterica.
Although most mutations were unique and appeared to be random in nature, a few of them occurred several times. For example, a mutation in the gene coding for flagellar transcriptional regulator flhD, resulting in a nonsynonymous amino acid change occurred in both serovars. FlhD protein is part of a transcriptional activator complex FlhD 4 C 2 , a master regulator of flagellar expression that then eventually allows the cell to be motile or non-motile based on the external and/ or internal stimuli 37 . In addition to FlhD, mutations in gene coding for flagellar motor switch protein (FliG) and flagellar basal body rod protein (FlgG) were also noticed. Physico-chemical environment, including temperature has been shown to affect flagellar expression and in general, the expression level is reduced with increased temperature. For example, Walker et al. 38 have shown reduced flagellar expression in S. Enteritidis grown at 37 °C in comparison to 20 °C. Similarly, Sirsat et al. 39 have shown the downregulation of flagellar genes of S. Typhimurium when exposed to sub-lethal thermal stress conditions through transcriptional profiling. In light of these reports and considering the fact that flagellar expression is a complex phenomenon that is intrinsically controlled at multiple levels 37 , we speculate that the genetic changes noticed in the genes related to flagellar expression would have resulted in its down regulation. In addition, mutations resulting in amino acid changes in other proteins such as NADP(H)-dependent aldo keto reductase, aspartate-semialdehyde dehydrogenase, CRISPR-associated helicase/ endonuclease Cas3 and glycosyl transferase (gtr) family I protein were also observed. Horizontally acquired gtr has been shown to contribute to the O-chain modification in the lipopolysaccharide component of outer membrane to enhance survival under adverse conditions in a host 40 . Thus, it is plausible that under exposure to stress conditions, the mutation induced in this gene would have resulted in the alteration of the composition of the O-chain to aid the survival of the organism. Additional research would be needed to confirm the precise role of the mutations observed in the genes and/ or intergenic regions.

Conclusion
Exposure of S. Agona and S. Mbandaka to repeated heat treatments (90 °C for 5 min) in a low water activity and high fat matrix did not result in increased fitness of the strains. However genetic changes were introduced which resulted in a population of isolates with a maximum of 28 allelic differences (38 SNPs) for S. Agona and 8 allelic (19 SNPs) differences for S. Mbandaka. This phenomenon could also potentially contribute to the higher genetic differences that can be observed between isolates sampled over a rather short timeframe and this knowledge is key during the interpretation of WGS results in a source tracking investigation involving isolates that could have been exposed to heat and/ or dry stress conditions.