A collection of 166 Listeria monocytogenes isolates from Canada and Switzerland were used in this collaborative study. This included: (i) 159 food and food processing environment isolates from Canada (n = 139) and Switzerland (n = 20), (ii) six isolates from sporadic human listeriosis cases in Switzerland, and (iii) one isolate from an asymptomatic human (Table S1). All human isolates were anonymized and no ethical approval was required as per the institutional and national guidelines. Isolates were stored at −70°C in brain heart infusion broth (BHIB, Difco, Fisher Scientific, Canada) +20% glycerol and routinely cultured at 30°C on BHI agar (BHIA, Difco, Fisher Scientific) plates.

Genomic DNA was isolated using the PureLink Mini Kit from Life Technologies, Canada. PicoGreen quantification was performed (Invitrogen, Canada) and DNA was assessed using the NanoDrop 2000 (Fisher Scientific). Genomic DNA samples of sufficient quality and quantity were sequenced by Genome Quebec (Montréal, QC, Canada) using TruSeq automated library preparation (Illumina) and paired-end, 100 bp sequencing on the Illumina Hi-Seq. Between 4.9 and 16.5 million high quality reads remained after quality control for each genomic library. Raw FASTQ files were trimmed using Cutadapt in Trim Galore! version 0.4.1 and de novo genome assembly was performed using SPAdes version 3.1.0 (careful option used; Bankevich et al., 2012). Low coverage (<10) and small contigs (<200 bp) were removed from assemblies using a custom perl script. Assemblies were subsequently annotated using Prokka version 1.5.2 (genus Listeria, species monocytogenes; Seemann, 2014). Assembled sequences were deposited into the NCBI Whole Genome Shotgun (WGS) database under Bioproject PRJNA329415.

To classify isolates into genetic lineages, a reference free, k-mer based single nucleotide variants (SNV) phylogeny was generated using the kSNP 3.0 program (Gardner et al., 2015) and reference isolates for the major lineages of L. monocytogenes (LI—F2365; LII—EGD-e; LIII—HCC23; LIV—J1-208). The resulting maximum parsimony tree (based on the consensus of 100 trees) clearly segregated the four lineages.

To group isolates based on their epidemiological context, in silico MLST was performed using the Center for Genomic Epidemiology's MLST typing tool (https://cge.cbs.dtu.dk/services/MLST/). Clonal Complexes (CCs) were assigned based on the Pasteur Institute schema (http://www.pasteur.fr/recherche/genopole/PF8/mlst/Lmono.html). Novel sequence types (STs) were confirmed using Sanger sequencing and submitted to the Pasteur Institute Database for new assignments (http://bigsdb.pasteur.fr/listeria/listeria.html).

Antibody-based serotyping was conducted on a subset of isolates (n = 91) within both the current study and previous studies (Arguedas-Villa et al., 2010; Kovacevic et al., 2013). Remaining isolates were assigned one or more possible serotypes by performing a MegaBLAST search (>95% nt identity) ncbi-blast+ v. 2.3.0 available at: ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/) for four genes used in a multiplex PCR developed by Doumith et al. (2004). Additionally, predictive serotypes were assigned to isolates with STs that are known to be associated with a specific serotype.

The genes and genomic regions evaluated in this study were (1) the plasmid replicon gene repA, used to indicate the presence of a plasmid (Kuenne et al., 2010), (2) emrE, representing Listeria genomic island 1 (LGI1), a 50 kb island with putative roles in stress tolerance and persistence (Gilmour et al., 2010), and (3) stress survival islet 1 (SSI-1), a five gene cluster previously identified as having a role in L. monocytogenes' response to cold, osmotic, and acid stress conditions (Ryan et al., 2010). Additionally, the coding sequence of inlA was investigated to determine if isolates possessed a full length sequence or a premature stop codon (PMSC) mutation. emrE, SSI-1, and inlA were screened for among isolate sequence assemblies using MegaBLAST (>95% nt identity) and repA was screened for using BLASTP (>30% aa identity over >80% coverage). inlA and repA sequences were then extracted from the isolate assemblies for further analysis.

Detection of repA sequences meant that at least one contig belonged to a putative plasmid. To identify additional plasmid associated contigs, isolate assemblies were aligned to the closed genome of L. monocytogenes EDG-e (Accession: NC_003210.1) using Contig Mover in Mauve version 2.3.1 (Rissman et al., 2009). Contigs not aligning to the EDG-e chromosome were compared to published L. monocytogenes plasmids (Kuenne et al., 2010) by BLAST. Contigs were excluded if they displayed open readings frames associated with chromosomal DNA (e.g., rRNA, tRNA) or did not align to any of the Listeria associated plasmids annotated by Kuenne et al. (2010): pLM33, pLM1-2bUG1, pLM5578, pLM80, and pLI100. A summary of the putative plasmid contigs found within each isolate can be found in Table S1.

Overnight cultures grown in BHIB at 37°C were standardized to 10⁹ CFU/ml using spectrophotometric methods, and diluted in pre-chilled BHIB to yield a final density of 10³ CFU/ml and stored at 4°C (previously described in Arguedas-Villa et al., 2010). The bacterial density was enumerated daily for the first 4 days and then bi-weekly for up to 5 weeks by plating on tryptic soy agar (BD, Fisher Scientific) +6% yeast extract (BD, Fisher Scientific). The resulting growth curves were fitted using a four parameter logistic model described by Dalgaard and Koutsoumanis (2001).

Isolates were assessed for salt and acid tolerance using modified versions of published protocols (Cotter et al., 2005; Van Der Veen et al., 2008; Bergholz et al., 2010). In short, overnight cultures grown in BHIB at 30°C were diluted in either BHIB+6% (w/w) NaCl or BHIB adjusted to pH 5 (with 1 M HCl) to achieve a final concentration of 10⁷ CFU/ml. From these cultures, 200 μl was added in duplicate (technical replicates) to 96-well-plates (Costar™ clear polystyrene, Fisher Scientific) that were incubated at 25°C in a microplate reader (Spectramax, V6.3; Molecular Devices, Sunnyvale, CA). A temperature of 25°C was used to assess isolate salt and acid tolerance under non-intracellular or cold stress conditions. The absorbance (A_600nm) of each well was recorded every 30 min until all isolates reached stationary phase (~26 h) and the resulting growth curves were fitted to the Baranyi and Roberts model (Baranyi and Roberts, 1994) using DMfit (v3.5) available on the ComBase browser (http://browser.combase.cc/DMFit.aspx).

Cultures grown for 24 h in BHIB at 20°C were diluted to 10⁷ CFU/ml in buffered peptone water (BD, Fisher Scientific) and 10 μl (10⁵ CFU) was spotted in duplicate (technical replicates) on the bottom of wells in lid-less 96-well-plates. The plates were then stored for 3 days at 20°C in desiccators (SICCO, Bohlender, Germany) pre-conditioned to 33% relative humidity (RH) using a saturated solution of MgCl₂ (protocol adapted from Hingston et al., 2015). The RH of the chambers was monitored throughout the desiccation periods using data loggers (included with desiccators). A temperature of 20°C was used to simulate desiccation conditions that might occur in a food plant. Following desiccation, the plates were rehydrated with 200 μl of BHIB, and incubated at 25°C in a plate reader where the A_600nm of each well was recorded every 30 min until all isolates reached stationary phase (~24 h). The resulting growth curves were then fitted to the Baranyi and Roberts model and the model parameters recorded.

For all four stress exposure experiments, a minimum of two biological replicates with two technical replicates each, were conducted for all isolates. Based on the findings of Aryani et al. (2015), the data was standardized for biological variability between replicates by dividing isolate growth parameters by the median value for each experimental run, thereby making the median equal to 1. The median was selected for standardization rather than the mean to avoid the influence of very stress sensitive isolates. Model parameters (LPD, lag phase duration; μ_max, maximum growth rate; N_max, maximum cell density) were averaged across biological replicates and presented as standardized (std) values. For isolates where the average std values had a standard deviation (SD) >0.05, additional replicates were completed to obtain more representative means. Isolates were considered tolerant or sensitive to cold, salt, or acid stress if they had an average std-μ_max > or < than 1 SD from the median, respectively. All remaining isolates were considered to have intermediate stress tolerance. For desiccation stress survival the model parameter of most interest was the LPD, indicating the time to (detectable) regrowth (TRG) that is negatively correlated with the number of cells, which survived the desiccation treatment. Isolates were classified as desiccation tolerant or sensitive if they had an average std-TRG < or > than 1 SD from the median, respectively. A standard curve generated using five cell levels (10¹–10⁵ CFU) produced a correlation of y = −0.25x + 2.07 (R² = 0.97) where y is the TDR and x is the log₁₀ number of viable cells in each well following desiccation.

To elucidate potential associations between the factors we investigated, statistical tests were performed using IBM SPSS Statistics version 23. Specifically, individual two-tailed T-tests and one-way ANOVAs with Tukey post-hoc tests were used to compare the average standardized stress tolerance model parameters of two (±plasmid, ±SSI-1, ±LGI1, lineage I and II, repA group 1 vs. group 2 isolates) or more groups (serotypes, CCs, inlA profiles, and sensitive, intermediate and stress tolerant groups), respectively. G^*Power 3.1.9.2 (Faul et al., 2007, 2009) was used to determine the minimum sample sizes required to ensure a power of 0.80 for all statistical tests. All data sets were accessed for outliers, homogeneity of variances (Levene's test), and normality (Shapiro-Wilk's test). Where homogeneity of variances could not be achieved, Welch T-tests and Welch ANOVAs in combination with Games-Howell post-hoc tests were used. P-values below 0.05 were considered significant for all comparisons.

Parsnp, a tool within Harvest suite of tools (Treangen et al., 2014), was used to perform core genome alignment of all 166 de novo assembled genomes and the reference L. monocytogenes EGD-e strain in order to identify single nucleotide variants (SNVs) within the core genome. SNVs clustered within 20 base pairs were removed as these may indicate repetitive regions containing more erroneous SNV calls. The remaining high quality SNVs were used to generate maximum likelihood trees using the RaxML version 8 (Stamatakis, 2014) on the CIPRES science gateway (Miller et al., 2010) using default parameters (including the GTRCAT nucleotide model and 100 bootstrap replicates). Corresponding heatmaps containing additional genotype and phenotype information were generated in R version 2.15.1 (Team, 2016) using the heatmap.2 function from the gplots library.

SNVs were also detected against the Listeria monocytogenes EGD-e (NC_003210.1) reference genome. SMALT version 0.7.6 (http://www.sanger.ac.uk/science/tools/smalt-0) with default parameters except “–i 330” was used to first align raw reads against the reference. Samtools version 1.2 (Li, 2011) was used on these assemblies to sort the aligned reads (“samtools sort”), remove potential PCR duplicates (“samtools rmdup”) and call the SNVs (“samtools mpileup”). Additional filtering of SNV calls included removing those with a read depth <50 and heterozygous genotypes (since our genomes are haploid) using the “bcftools filter” command. SNVs found in repetitive regions of the genome as assessed by the index of repetitiveness (Schwender et al., 2004) were also removed manually. The remaining high confidence SNVs were annotated using SNPEff version 4.1 (Cingolani et al., 2012) with the Listeria_monocytogenes_EGD_e_uid61583 annotation. Synonymous SNVs were also removed in the end for identification of non-synonymous or potential regulatory SNVs that may be contributing to phenotypic differences in cold growth.

SNPSift version 4.1 “CaseControl” (Cingolani et al., 2012) was used to run a Fisher Exact test to identify SNVs that were significantly associated with case vs. control groups. To identify SNVs only found in tolerant isolates, these were used as the case group, while all others were used as the control group. Since this did not yield any results, subsequently, the sensitive isolates alone were used as the control group so as to allow SNVs to be seen in intermediate growers. This method has certain limitations in that certain associations may require very large sample sizes to become statistically significant, especially considering genetic heterogeneity leads to the same phenotype or the potential for multiple SNVs to interact. An alternative approach, Random Forests™ (Breiman, 2001), was also used to discover important SNVs in distinguishing stress tolerant and sensitive groups since previous genome wide association studies have shown random forests outperform the Fisher Exact test in these special cases (Lunetta et al., 2004; Schwender et al., 2004; Bureau et al., 2005; De Lobel et al., 2010; Bulinski et al., 2011). The RandomForest™ version 4.6–10 library was used in R with default parameters except “importance = TRUE, proximity = TRUE, ntree = 5000.” This allows the method to run as a classifier that then ranks SNVs on their ability to classify isolates based on their phenotypic designation.

Annotated draft genomes were submitted to IslandViewer 3 (Dhillon et al., 2015) using L. monocytogenes EGD-e (NC_003210.1) as a reference for contig reordering. Genomic islands were predicted using IslandPath-DIMOB (Hsiao et al., 2005) and SIGI-HMM (Waack et al., 2006). Predicted genomic islands positioned on the genome within <10 kb of each other were merged into one single region.

To form groups of similar genomic islands, the genetic distance between genomic island sequences was computing using Mash (parameter –s 2000; Ondov et al., 2016) and groups of similar sequences were identified using hclust and cutree in R.

The complete sequenced genome assembly sizes of the isolates ranged from 2.56 to 3.13 Mbp with a mean size of 2.97 Mbp (Table S1). Isolates belonged to one of three different lineages: LI (n = 44, serotypes 4b, 1/2b, 3b, and 3c), LII (n = 121, serotypes 1/2a, 1/2c, and 3a), and LIII (n = 1, serotype 4c). The majority of isolates were serotype 1/2a (n = 92), followed by 1/2c and 4b (n = 25 each), 1/2b (n = 18), 3a (n = 2), and 3b and 4c (n = 1 each; Table 1). The exact serotype was not determined for two remaining isolates. Beyond serotypes, our isolates belonged to 36 different known STs and a further nine were assigned novel STs (ST1017-1025). Isolates also belonged to one of 29 different CCs with a further seven isolates being unique non-clonal singletons. The most prevalent CCs in the collection were CCs 9, 8, and 7 (Table 2). Other less common CCs in decreasing prevalence included CCs 11, 155, 1, 3, and 321 (Table 2). Interestingly, only one CC121 isolate existed in our collection. This is surprising given that CC121 is often highly prevalent among L. monocytogenes food-associated isolates (Parisi et al., 2010; Chenal-Francisque et al., 2011; Martín et al., 2014; Ebner et al., 2015; Maury et al., 2016).

The plasmid replication gene, repA, was observed in 55% (n = 92) of our isolates, with a prevalence of 41 and 61% among LI and LII isolates, respectively. Notably, one isolate was observed to contain two putative plasmids as indicated by the presence of two different repA containing contigs of 61 and 69 kb. Among serotypes, repA was present in 100% of 3a isolates, 84% of 1/2c isolates, 78% of 1/2b isolates, 53% of 1/2a isolates, and 16% of 4b isolates (Table 1). Among CCs, plasmids were observed in >80% of CC 3, 5, 9, 11, and 321 isolates (Table 2).

A phylogeny, constructed on repA sequences as described in Kuenne et al. (2010), divided the sequences into two groups. Group 1 included isolates from serotypes 1/2a, 1/2b, 1/2c, and 4b with estimated plasmid sizes ranging from 26 to 88 kb. Group 2 included 1/2a, 1/2b, 1/2c, and 3a serotype isolates, harboring significantly larger plasmids (p < 0.0005, 55–100 kb) than those from group 1. These sizes are in line with those observed in Kuenne et al. (2010), supporting the assertion that these contigs belong to plasmids. The most prevalent plasmid size (56553–56554 bp) was observed for 26 isolates from seven different CCs and from both lineages I and II. Also noteworthy is that isolates from Switzerland and Canada contained plasmids with 100% nucleotide identity.

Premature stop codons (PMSCs) in inlA were observed in 20% of our isolates encompassing seven (Table S1) of 19 published PMSCs (Jonquieres et al., 1998; Olier et al., 2002, 2003; Rousseaux et al., 2004; Nightingale et al., 2005, 2008; Orsi et al., 2007; Van Stelten and Nightingale, 2008; Van Stelten et al., 2010; Wu et al., 2016) and one novel PMSC at 760aa, which was identified in two serotype 1/2c isolates. The most common PMSC occurred at 9aa (n = 13) and was associated with CC9, serotype 1/2c isolates. Ten of the 14 remaining CC9 isolates also had one of four inlA PMSCs (326, 576, 685, 760aa) and all CC321 isolates contained inlA PMSC's at 700aa (Table S1). An additional 13 isolates, all from the serotype 4b CCs 6 and 315, contained a three codon deletion mutation previously reported in Kovacevic et al. (2013). With the exception of CCs 5 and 9, all isolates from the same CC either contained full length inlA or a truncated version.

During the screening of the whole genome sequences, the absence of lmo1078 was noted among serotype 4b isolates. This gene, which encodes a UDP-glucose pyrophosphorylase, has been previously demonstrated to have a role in L. monocytogenes cold growth (Chassaing and Auvray, 2007). It was also observed that 70% (n = 116) of strains possessed SSI-1 with this island being most prevalent among serotype 1/2c isolates (100%) followed by 1/2b (94%), 1/2a (71%), and 4b (16%; Table 2). Furthermore, all isolates from CCs 3, 5, 7, 8, 9, 155, 224, 315, and 321 contained SSI-1 (Table 2). All remaining isolates possessed a homolog to F2365_0481 in place of SSI-1 (Ryan et al., 2010), with the exception of the CC121 isolate which possessed lin0464 and lin0465 homologs as reported in Hein et al. (2011).

The LGI1 indicator gene, emrE, was found in 16 of our isolates and as previously reported (Gilmour et al., 2010; Althaus et al., 2014; Kovacevic et al., 2015), all originated from Canada and 14 were serotype 1/2a ST120-CC8. The remaining two isolates represented novel STs (ST1022 and 1025) that also belonged to CC8. All emrE containing isolates also harbored SSI-1 and full length inlA.

All L. monocytogenes isolates were evaluated on their ability to grow in cold (4°C), salt (6% NaCl), and acid (pH 5) stress conditions as well as survive desiccation stress (33% RH). The cold growth plate count data was modeled using the Dalgaard and Koutsoumanis (2001) logistic model because it was more accommodating of fewer sampling points [average R² = 0.998, mean standard error (MSE) = 0.129]. From the std-μ_max values, 13 isolates were classified as cold sensitive and 18 were classified as cold tolerant with average std-μ_max values of 0.85 ± 0.08 and 1.09 ± 0.02, respectively (Figure 1A). For the salt, acid, and desiccation tolerance assays, the Baranyi and Roberts (1994) was suitable for modeling the spectrophotometrically obtained data with average R² and MSE-values ranging from 0.997 to 0.998 and 0.003 to 0.017, respectively. Overall, 27 and 17 isolates were classified as salt sensitive and tolerant with average std-μ_max values of 0.83 ± 0.05 and 1.16 ± 0.05 (Figure 1B); 26 and 22 isolates were classified as acid sensitive and tolerant with average std-μ_max values of 0.64 ± 0.14 and 1.34 ± 0.12 (Figure 1C); and 20 and 23 isolates were identified as desiccation sensitive and tolerant isolates with average std-TRGs of 0.81 ± 0.06 and 1.22 ± 0.11 (Figure 1D), respectively.

Five isolates were classified as sensitive to three out of four stresses and 10 isolates were sensitive to two stresses, seven of which were salt and acid sensitive (Figure 2A). Only two isolates were classified as tolerant to three out of the four stresses and another 14 isolates were tolerant to two of the four stresses (Figure 2B). Twenty additional isolates displayed a total of 16 combinations of overlapping sensitive and tolerant phenotypes (Table S1). The most common overlapping phenotypes were salt and acid sensitive (n = 10), salt sensitive and desiccation tolerant (n = 6), cold and acid tolerant (n = 5), and cold tolerant and salt sensitive (n = 5).

Isolates designated sensitive, intermediate or tolerant to one stress were analyzed to determine if they significantly differed in their tolerances to other stresses. When grown in 6% NaCl, acid tolerant isolates had larger std-μ_max values compared to acid sensitive isolates (p = 0.03, x 1.02 vs. 0.95). Similarly, salt sensitive isolates had smaller std-μ_max values (x = 0.80) in BHIB pH 5 compared to intermediate (p < 0.0005, x = 1.03) and salt tolerant isolates (p = 0.006, x = 1.00).

Between lineages, the only significant difference observed was that LI isolates had significantly larger std-μ_max values in BHIB pH 5 than LII isolates (p < 0.0005, x = 1.13 vs. x = 0.94). Additional significant differences were observed between serotypes. At 4°C, serotype 1/2a isolates had significantly larger (p = 0.017) std-μ_max values compared to serotype 1/2c isolates (Figure 3). In support of this, serotype 1/2a isolates accounted for 61% of the cold tolerant isolates and only 31% of cold sensitive isolates compared to a 55% prevalence of this serotype in the collection (Table 1). Similarly, serotype 1/2c isolates accounted for 38% of cold sensitive isolates despite a 15% overall prevalence in the collection (Table 1). When isolates were grown in 6% NaCl, no significant differences were observed between serotypes (Figure 3), however, 71% of salt tolerant and 67% of salt sensitive isolates were serotype 1/2a isolates, relative again to a prevalence of 55% in the collection (Table 1).

In BHIB pH 5, serotype 1/2a isolates had significantly smaller (p = 0.027) std-μ_max values than serotypes 1/2b, 1/2c, and 4b (Figure 3). In agreement with these findings, 85% of acid sensitive isolates were serotype 1/2a whereas only 32% of acid tolerant isolates were serotype 1/2a (Table 1). No significant differences (p > 0.05) were observed between serotypes with respect to desiccation stress std-TRGs (Figure 3).

Beyond the stress tolerances of lineages and serotypes, some significant differences were also detected between CCs. As a minimum of six isolates per CC were needed to ensure a power >0.80 for ANOVA results, statistical analyses were only performed using CCs 1, 3, 5, 6, 7, 8, 9, 11, 155, and 321. Figure 4 shows the average levels of cold (std-μ_max), salt (std-μ_max), acid (std-μ_max), and desiccation (std-TRG) tolerance among CCs with three or more isolates. At At 4°C, no significant differences were found between the growth rates of different CCs, however, it was interestingly to see that CCs associated with 4b isolates had both the lowest and highest average average std-μ_max values at 4°C, demonstrating why stress tolerance differences were not observed between this serotype and others at 4°C (Figure 4A).

In 6% NaCl, CC7 (1/2a) isolates had significantly (p < 0.05) smaller std-μ_max values compared to CCs 5 (1/2b), 8 (1/2a), 11 (1/2a), and 155 (1/2a; Figure 4B). This highlights the range of salt tolerances between CCs within the same serotype and again explains why no significant differences were observed at the serotype level for salt tolerance. In support of the results shown in Figure 4B, 67% of CC2 isolates were salt tolerant while 50% of CC224, 47% of CC7, and 22% of CC9 isolates were salt sensitive (Table 2).

In BHIB pH 5, CC5 (1/2b) isolates exhibited significantly (p < 0.05) larger std-μ_max values than CCs 7, 155, and 321, and CC321 isolates additionally had smaller (p < 0.05) std-μ_max values compared to CCs 1, 3, and 11 (Figure 4C). CC1 (4b) isolates also had significantly (p = 0.02) smaller std-μ_max values compared to CC7 (1/2a) isolates. From Figure 4C it can be seen that lineage I isolates were more acid tolerant than LII isolates, as the five CCs with the highest average std-μ_max values in BHIB pH 5 were all from LI while the five CCs with the lowest average std-μ_max values belonged to LII (predominantly 1/2a isolates). Notably, 57% of CC5 and 67% of CC4 isolates were acid tolerant while 83% of CC321, 67% of CC20, and 29% of CC7 isolates were acid sensitive (Table 2).

No significant differences were found between the desiccation stress std-TRGs of different CCs (Figure 4D). Nevertheless, CC224 (1/2b) had the smallest average std-TRGs and correspondingly, 75% of these isolates were classified as desiccation tolerant. CC11 (1/2a) had the next smallest std-TRGs while CCs 1 and 4 (both 4b) had the two largest average std-TRGs.

Although plasmids were identified in 55% of all isolates, a higher percentage of plasmid carriers were observed among acid tolerant (73%), desiccation sensitive (75%), and desiccation tolerant (60%) isolates as compared to cold tolerant (33%) and acid sensitive (46%) isolates (Figure 5). Within LII, plasmid-positive isolates had smaller std-μ_max values at 4°C (p = 0.024, x = 1.00 vs. 1.02) and larger std-μ_max (p < 0.0005, x = 1.01 vs. 0.86) values when grown in BHIB pH 5 compared to their plasmid-free counterparts. No significant differences were found between the stress tolerance levels of LI plasmid-harboring and plasmid-free isolates. It was, however, observed that isolates containing repA group 1 plasmids; which were significantly (p < 0.0005) smaller than group 2 plasmids, had smaller std-μ_max (p = 0.002, x = 0.98 vs. 1.04) values in 6% NaCl.

Full length inlA was observed in 72% of isolates, where a higher percentage of the intact gene prevailed among cold (89%), salt (94%), and acid (82%) tolerant isolates while a lower prevalence was detected among desiccation tolerant isolates (35%; Figure 5). Statistically, isolates with full length inlA had significantly larger std-μ_max values at 4°C than isolates with an inlA PMSC (p = 0.001, x = 1.01 vs. 0.97). Additionally, serotype 4b isolates possessing a three-codon deletion in inlA, had significantly shorter desiccation stress std-TRGs (p = 0.002 x = 0.94 vs. 1.05) compared to serotype 4b isolates with full length inlA. No significant associations were found between inlA profiles and salt or acid stress tolerance.

In 6% NaCl, isolates containing SSI-1 had significantly smaller std-μ_max values (p = 0.004, x = 0.98 vs. 1.03) than isolates without SSI-1 though this difference was not large. To determine potential associations between LGI1 harborage and stress tolerance, serotype 1/2a isolates containing the island were compared to other 1/2a LGI1-negative isolates but similar to SSI-1, no significant differences in stress tolerances were detected.

Figure 6 shows a whole genome SNV phylogeny of all 166 L. monocytogenes isolates with their corresponding genetic and phenotypic properties. In this figure, groups of closely related isolates that share the same phenotypes can be seen. However, also shown are several cases where neighboring isolates have opposing stress tolerances. Of particular interest was whether specific SNVs could be related to isolates possessing the same stress tolerance phenotypes, however, none were detected to be uniquely shared among stress tolerant isolates that weren't also seen in intermediate or sensitive isolates. Among stress sensitive isolates, unique SNVs shared by subsets of isolates were identified, but no single SNV was prevalent among >4 isolates from the same stress sensitive phenotype group. In contrast, a large number of SNVs were uniquely observed for one or two isolates from the same stress sensitive group, causing frameshifts, premature stop codons, loss of start codons, or missense variants. Information regarding the SNVs identified among all sensitive and tolerant isolates are presented in Tables S2–S9. Notably, a number of stress sensitive isolates contained different PMSCs in several σ^B regulator genes. A cold and desiccation sensitive isolate contained a PMSC in rsbS as did two other desiccation sensitive isolates. Furthermore, an additional cold sensitive isolate contained a PMSC in rsbV, and two desiccation sensitive isolates contained PMSCs in rsbU.

All L. monocytogenes isolates were predicted to harbor 1,318 genomic islands in total, resulting in an average of eight genomic islands per genome. These islands were clustered into 200 groups of similar sequences. The conservation of genomic island groups across L. monocytogenes lineages ranged from unique to a single isolate to conserved in 97 isolates (58%). Most frequently, intermediate groups conserved in subsets of monophyletic isolates were observed, as can be expected from a combination of vertical inheritance of the genomic islands with further modifications by mutation, insertions, and deletions. The clustering of L. monocytogenes isolates based on the presence or absence of groups of genomic islands reflected their phylogenetic proximity and did not relate particular genomic island content in isolates with the exhibition of similar phenotypes. Indeed, no single genomic island was found to occur in a large proportion of strains with a given phenotype (Figure S1).

An important finding from this study was that closely related isolates from within the same clonal complexes exhibited opposing stress tolerances, suggesting that minor genetic differences can also exert great impact on stress tolerance phenotypes. This was observed in a study by Hoffmann et al. (2013), where a single thymine deletion in the σ^A-like promoter region of betL, encoding an osmolyte transporter specific for betaine uptake, dramatically increased betL transcription, and hence the osmo- and chill-tolerance. Karatzas et al. (2003) reported that a spontaneous high hydrostatic pressure tolerant L. monocytogenes mutant of Scott A, contained a single codon deletion in ctsR, a negative regulator of several heat shock and general stress proteins, that also conferred increased thermo-tolerance and resistance to H₂O₂. Additionally, in Metselaar et al. (2015) a number of spontaneous acid tolerant mutants were found to contain SNVs in the ribosomal protein gene rpsU. None of these aforementioned mutations were detected in the present study.

Analyses to determine if unique SNVs could be detected among isolates from individual stress tolerance phenotype groups were also performed. Studies identifying the genetic basis of phenotypic traits using the variation within natural populations are known as a genome-wide association studies (GWAS). While GWAS have been effective for identifying mutations responsible for phenotypic traits in humans, the clonal nature of bacterial replication where mutations can reach a high frequency on a single genetic background, makes it difficult to distinguish mutations responsible for an observed phenotype (Read and Massey, 2014; Falush, 2016). As a result, bacterial molecular epidemiology has focused on identifying clonal lineages with particular phenotypic properties rather than identifying the specific genetic variants responsible. Recently, Earle et al. (2016) used a GWAS approach to successfully identify genes and genetic variants underlying resistance to 17 antimicrobials in over 3000 isolates of taxonomically diverse clonal and recombining bacteria. While these results show the potential of bacterial GWAS, antimicrobial resistance is usually gained during antimicrobial exposure and thus it is more likely that the traits evolve on multiple independent backgrounds making them easier to detect (Falush, 2016). To date, the identification of mutations responsible for more complex phenotypes such as those evaluated in the present study, remain challenging.

Our analyses did not detect any unique SNVs among more than one isolate from the same stress tolerant group, suggesting homoplasy among stress tolerant phenotypes where mutations evolve independently to confer tolerance. On the contrary, a few different SNVs were identified among four or fewer isolates from the same stress sensitive groups. Notably, the cold sensitive phenotypes of two isolates may be associated with PMSCs detected in the σ^B regulator genes rsbS and rsbV (Voelker et al., 1995). In support of this hypothesis, deletion of rsbV in a L. monocytogenes mutagenesis study resulted in impaired cold stress tolerance (Chan et al., 2008). Three desiccation sensitive isolates also contained different PMSCs in rsbS, including one isolate which was also cold sensitive. An additional two desiccation sensitive isolates contained different PMSCs in another σ^B posttranscriptional regulator, rsbU. Given the importance of σ^B in L. monocytogenes' adaptation to several environmental stresses (Wiedmann et al., 1998; Kazmierczak et al., 2003), it is possible that these mutations contributed to the reduced desiccation tolerance of these isolates. In fact, in Huang et al. (2015) an L. monocytogenes sigB mutant demonstrated reduced desiccation survival relative to the wildtype strain.

In this study, specific genomic islands could not be exclusively associated with a particular stress tolerance phenotype. This is to be expected as genomic islands often contain virulence factors, and in general these are overrepresented in genomic islands as compared to the chromosome (Sui et al., 2009). In L. monocytogenes in particular, genomic islands have been associated with virulence, heavy metal resistance, and benzalkonium chloride efflux (Gilmour et al., 2010; Kuenne et al., 2010; Kovacevic et al., 2015). SSI-1 as described above, has been previously associated with L. monocytogenes' stress response, but was not significantly correlated with the stress tolerance phenotypes examined in the present study.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.