Complete genome of the barotolerant Listeria monocytogenes RO15 strain and comparison with other strains isolated from food and food processing environments

Background: High pressure processing (HPP; i.e. 100 - 600 MPa pressure depending on product) is a non-thermal preservation technique adopted by the food industry to decrease significantly foodborne pathogens, including Listeria monocytogenes , from food . However, susceptibility towards pressure differs among diverse strains of L. monocytogenes and it is unclear if this is related to their genomic content. Here, we tested the barotolerance of 10 different L. monocytogenes strains, from food and food processing environments and widely used reference type strains, to pressure treatments with 400 and 600 MPa. Genome sequencing and genome comparison of the tested L. monocytogenes strains were performed to investigate the relation between genomic profile and pressure tolerance. Results : None of the tested strains were tolerant to 600 MPa. A reduction of more than 5 log 10 was observed for all strains after 1 minute 600 MPa pressure treatment. L. monocytogenes strain RO15 showed no significant reduction in viable cell counts after 400 MPa for 1 minute and was therefore defined as barotolerant. Genome analysis of so far unsequenced L. monocytogenes strain RO15, 2HF33, MB5, AB199, AB120, C7, and RO4 allowed us to compare the gene content of all strains tested. This revealed that the three most pressure tolerant strains had more than one CRISPR system with self-targeting spacers. Furthermore, several anti-CRISPR genes were detected in these strains. Pan-genome analysis showed that 10 prophage genes were significantly associated with the three most barotolerant strains. Conclusions: L. monocytogenes strain RO15 was the most pressure tolerant among the selected strains. Genome comparison suggests that there might be a relationship between prophages and pressure tolerance in L. monocytogenes . Analysis of variance (ANOVA), average nucleotide identity (ANI), coding sequence (CDS), colony forming units (cfu), high pressure processing (HPP), megapascals (MPa), multilocus sequence typing (MLST)

enterica [18], and L. monocytogenes [19,20]. In L. monocytogenes strains, pressure tolerance varies between 300 -500 MPa [19]. To our knowledge, genomic profiling and comparison has not been performed for barotolerant and barosensitive L. monocytogenes strains. In this study, we selected 10 L. monocytogenes strains, isolated from either foods, food processing environments or clinical sources, to test their tolerance towards HPP at 400 and 600 MPa, compared their genomes, and investigated whether genetic traits may be associated with pressure tolerance.

Pressure treatment and reduction of viable cell counts
Reduction of viable cell counts (colony forming units, cfu) after pressure treatment at 400 and 600 MPa for 1 minute showed that pressure tolerance differs between strains. The variance in colony forming units (cfu) of the treated samples was significantly larger (p<0.02, see ANOVA results below) compared to control, indicating a high variance in the level of tolerance of L. monocytogenes. At 400 MPa for 1 minute, the 10 strains exhibited an average log 10 reduction of 0.57 cfu/ml, ranging from 0.05 log 10 cfu/ml for the most barotolerant strain (RO15) to 2.07 log 10 cfu/ml for the most pressure sensitive strain (EGD-e). Similar results were obtained at 600 MPa. Here, the 10 strains exhibited an average log 10 cell number reduction of 7.06 cfu/ml, with a range of 5.42 to 8.27 log 10 cfu/ml (Table 1). Strain RO15 was also the most barotolerant strain based on an initial screening including several other L. monocytogenes strains (Supplementary text 1, Table S1, Figure S1).
One-way analysis of variance (ANOVA), comparing the pressure (400 MPa) treated samples of each strain to the control samples, showed that the mean cfu/ml for treated samples was significantly lower than the control for all strains (p<0.02), except for strain RO15 (p=0.15) ( Figure 1a). The same statistical analysis for HPP at 600 MPa indicated that all strains including strain RO15 had significantly (p<0.01) lower cfu/ml in treated samples compared to controls (Figure 1b).

Genome sequencing, general features and RNA-Seq
We sequenced the genomes of seven L. monocytogenes strains for which no genome sequences were available in public databases. Genomes of strains RO15, 2HF33, MB5, AB199, AB120, C7, and RO4 were sequenced using Illumina MiSeq instrument and strain RO15 was sequenced using Pacbio RSII. Sequences of these seven strains were assembled and three additional strains (ScottA, F2365 and EGD-e, genome sequences of which are available in public databases) were used for comparative genome analysis. Assembly of PacBio long reads of strain RO15 resulted in one continuous 3,042,507 bp sized chromosome ( Figure 2) at an average 308-fold sequencing coverage. In addition to the chromosome, a contig with a complete circular prophage sequence of 38,811 bp was obtained, which was also found as part of the chromosome (2729417 -2770759 bp). All genomes had a similar size, GC content and number of coding sequences (CDS) ( Table 2).
The serotype of the strains was reported in previous studies [21,22,23,24,25] (Table 2). Although RO15 was reported as a serovar 4b strain using multiplex PCR in a previous study [21], our genome-based prediction suggests that it belongs to serovar 1/2a. Sequences for ORF2110 and ORF2819 primers used for identification of serovar 4b strains [24] were not found in the genome of strain RO15 based on our sequence match analysis, whereas sequences for primers targeting lmo0737, which are indicators of serovar 1/2a strains [24] were present. Multilocus sequence typing (MLST) based on 7 loci, core genome MLST (cgMLST; based on Moura scheme [26]), and clonal complex (CC) of the strains were also assigned based on the genome sequences (Table  2). Lineage assignment showed that only ScottA and F2365 belonged to lineage I and all other strains were members of lineage II (Table 2).
In addition to genome sequencing, transcriptome analysis using RNA-seq for strains RO15 and ScottA provided a basic view of transcriptional activity of the genomes. Differential expression analysis indicated that virulence genes and heat shock genes were upregulated after high pressure treatment in strain ScottA (Table S2).

Methylated DNA Motifs
High PacBio sequencing coverage allowed us to analyze DNA methylation modification motifs in strain RO15 in addition to assembling the genome. Using the long read modification data, 11 methylated sequence motifs were detected in the genome ( Table 3). None of the detected motifs had a partner motif, i.e. a reverse-complementary sequence, but all detected motifs were only partially modified in the genome with less than 50% methylated motifs. While most of the motifs have already been deposited in the REBASE database [29] DADGYATYA, WNNTVVGCNTWNH, AHNBAACA, AGNNARNWW were novel, i.e. they have not been described as potential methylation sites previously. None of the detected motifs have been reported as recognition sequence motif for a restriction enzyme in the REBASE database.
We predicted one type II cytosine-5 DNA methyltransferase gene (OCPFDLNE_00657), and three type II N4cytosine or N6-adenine DNA methyltransferase genes (OCPFDLNE_02168, OCPFDLNE_02626, OCPFDLNE_02808) in the strain RO15 genome. In addition, genes for a type IV methyl-directed restriction enzyme (OCPFDLNE_00324) and type II restriction enzymes (OCPFDLNE_00658, OCPFDLNE_02625, OCPFDLNE_02807) were predicted. These predicted methyltransferases and restriction enzymes had significant (e-value <1E-50) BLASTP hits in the REBASE protein sequences database [29]. However, their recognition sequences are unknown. OCPFDLNE_02625 and OCPFDLNE_02807 type II restriction enzyme genes in prophage regions were expressed based on RNA-seq data.

Genomic comparison of the selected strains.
In the genome assembly we also identified a circular phage as an independent contig. While looking after phageoriginating genome parts, prophage prediction revealed five prophage regions (10.7 kb to 47.9 kb sized) in the RO15 genome ( Figure 2). Prophage (region) 5 had the same sequence as the separate circular phage. RNA-Seq count data suggested that most of the genes of prophage 1,3 and 5 in RO15 were transcribed. Subsequently, the prophage prediction was performed for all strains, which showed that all strains had at least one or more prophage regions in the genome with a maximum of six prophage regions in strain C7 ( Table 4). One of the predicted prophage regions in ScottA was also transcribed based on the RNA-seq data.
Genome alignment of the selected nine strains against the strain RO15 ( Figure 3) showed large sequence gaps between strain RO15 and the other strains. The gaps were mainly related to prophage regions with prophage region 3 being specific for strain RO15. Prophage region 2 was only seen in the strain 2HF33. Prophage region 4 was partially seen in strains 2HF33, EGD-e, MB5 and ScottA. Similarly, prophage region 5 was partially seen in strains 2HF33, C7, and RO4. In addition, strains ScottA and F2365 had less aligned regions compared to the other strains. Average nucleotide identity (ANI) results showed that ANI scores of strains ScottA and F2365 with other strains were lower than 95 (Table S3). To study how the large sequence gaps are localized in the other publicly available genomic sequences of L. monocytogenes, all 200 complete L. monocytogenes genomes in RefSeq database [30] were downloaded and aligned to strain RO15 ( Figure S2). The visualization of the alignment showed that prophages 2, 3 and 4 were seen only in the minority of the L. monocytogenes strains, while prophages 1 and 5 were seen in the majority of the L. monocytogenes strains. The comparison indicated that these prophage regions were common locations of variance between L. monocytogenes strains.
Annotation of CRISPR systems revealed that the RliB-CRISPR system was seen in all the strains, which is in line with previous RliB-CRISPR system studies [31,32]. Interestingly, CRISPR I or CRISPR II system genes were present only in strains RO15, 2HF33, and MB5 (Figure 4), which exhibited the lowest reduction in log 10 cfu/ml with the 400 MPa pressure treatment ( Table 1). Number of spacers in these strains were also significantly (p < 0.001) higher compared to barosensitive strains. The alignment of the spacer sequences of the CRISPR systems back to the genome itself revealed that only strains RO15, 2HF33, and MB5 contained self-targeting spacers with 100% identity. As expected, spacer sequences were aligned to the prophage regions, except one spacer sequence in RO15, which aligned to the addB gene (OCPFDLNE_02427) located in the chromosome encoding an ATP-dependent helicase.
As antibiotic resistance and pressure tolerance were linked in a previous study [19], we also searched antibiotic resistance genes in all strains to identify relation of antibiotic resistance genes and pressure tolerance. The same antibiotic resistance gene families were detected in all the 10 strains (Table S4). Multiple sequence alignment of the detected antibiotic resistance genes and the following average distance tree generation showed that amino acid sequence of the norB gene (encoding for a quinolone resistance protein) (OCPFDLNE_03068) was slightly different in barotolerant strain RO15 compared to the other strains ( Figure S3). Similarly, for the lin gene (encoding for a lincomycin resistance protein) (OCPFDLNE_00980) amino acid difference was seen for strains 2HF33 and RO15 ( Figure S3).
Roary pan-genome pipeline suggested that the 10 L. monocytogenes strains contain a total of 4820 orthologous gene clusters. Of these genes, 2247 were core genes found in all strains and 2573 genes were accessory genes found in at least one strain. The core genome was used for constructing a phylogenetic tree, which indicated that the two serovar 4b strains (ScottA and F2365) are closely related to each other and cluster separately from the serovar 1/2a strains in the phylogenetic tree. In addition, there was also a clear difference in accessory genome for serovar 4b strains compared to serovar 1/2a strains ( Figure 5).
To test any significant association between pressure tolerance and (clusters of) genes in the accessory genome, we performed a pan-genome wide association analysis using Scoary. Based on the pairwise comparison of reduction in log 10 (cfu/ml) against strain RO15 using Student's t-test (Table S5), a statistically significant difference compared to RO15 was seen in all strains except MB5 and 2HF33. Therefore, strains RO15, 2HF33, MB5 were used as barotolerant strains for the pan-genome association analysis. Scoary results showed that 13 gene clusters (Table S6) had p-values < 0.01 for the association with barotolerant strains, i.e., the genes were only seen in barotolerant strains. Of these, 10 gene clusters were located in the prophage regions, but most of the genes were annotated as hypothetical proteins (Table S6). These prophage genes were also searched for in 200 complete L. monocytogenes strains and found in some strains (Table S6). Interestingly, genes OCPFDLNE_02579 and OCPFDLNE_02580 were only seen in strains that harbored anti-CRISPR genes, Cas type I or II CRISPR system and self-targeting spacers (Table S6).

Discussion
HPP is commonly used in the food industry to inactivate foodborne pathogens and spoilage organisms. Depending on the type of food product, up to ~600 MPa pressure is applied [9]. It has also been shown that combinations of HPP and biocontrol agents, such as virulent bacteriophages and bacteriocins, have synergistic effects in the inactivation of the target pathogens in food [34,35]. The pressure tolerance level varies between different species, strains, and even isolates [16,19,36,37]. Here, we tested pressure tolerance of 10 L. monocytogenes strains, compared their genomes and predicted genome features related to pressure tolerance.
Pressure treatment using 600 MPa for 1 minute caused more than 5 log 10 reduction in all selected strains, which is generally considered sufficient for food safety regulations. According to the Food Safety and Inspection Service (United States Department of Agriculture), 5 log 10 reduction is considered a full lethality treatment for L. monocytogenes [38]. This suggests that the strains used in this study were relatively sensitive to 600 MPa pressure treatment. Similarly, an earlier study reported that pressure treatment at 500 MPa for 10 minutes provided sufficient reduction in viable counts to reach desired level of safety for all except one L. monocytogenes strains tested [19]. However, the lower pressure levels, such as 400 MPa, would be more relevant for industrial applications, therefore, we focused more on the results obtained with the 400 MPa treatment. Among the strains tested, RO15 was the most barotolerant strain when processed at both 400 and 600 MPa, which is in line with the initial selection process of strains (Supplementary text 1). Based on ANOVA, RO15 was the only strain for which no statistically significant reduction was observed with the 400 MPa treatment (Figure 1a). Therefore, we defined the strain RO15 as barotolerant. It has been shown that the level of inactivation by HPP can be affected by the chemical composition of the food product. Some food products, such as milk and cheese, have baroprotective effects on microorganisms during HPP treatment [39,40]. However, experiments in this study were performed using TSBYE broth and thus do not take into account the effects of food composition on HPP. Temperature during pressure treatment can have an effect on microbial inactivation depending on treatment time, and it has been shown that increasing the temperature from 25 °C to 35 °C has a significant effect on the viability loss of L. monocytogenes [41]. In this study, temperature during pressure treatment increased up to 33 °C for 400 MPa treatment and 38 °C for 600 MPa treatment due to adiabatic heating. However, since treatment time was only one minute, the effect of temperature rise was expected to be minimal.
Serotypes of the selected strains were shown in the earlier studies based on agglutination method, multiplex-PCR and genome sequencing [21,22,23,24,25] (Table 2). Here, we predicted the serotype of all strains based on their genome sequences. This confirmed the previous serotyping of all strains except RO15, which was reported to be a serotype 4b strain based on multiplex-PCR [21]. Serotype 4b strains should contain serotype 4 marker genes (ORF2110 and ORF2819) and lack lmo0737 marker gene sequence as described by Doumith et al. [24]. However, according to the genome sequence, serotype 4 marker genes were absent, and a lmo0737 homologue was present suggesting that the serotype of strain RO15 is indeed 1/2a. The phylogenetic tree that we created based on the core genome showed that there is considerable genetic variation between serovars. Serovar 4b strains (ScottA and F2365) clustered separately from the other strains ( Figure 4). Strain RO15 was not clustered together with other serovar 4b strains hence these results supported our gene-based prediction of serotype 1/2a of strain RO15.
PacBio sequencing provides not only the genome sequence but also methylation data [42], which gives an opportunity to analyze restriction-modification systems and their recognition motifs. In this study, we did not observe a genuine methylated motif based on PacBio sequencing data of strain RO15. Similarly, PacBio methylation data of nine other published L. monocytogenes strains without genuine methylated motifs were also seen in REBASE PacBio list database [29], which showed that it is not uncommon to have only partially methylated motifs in L. monocytogenes.
Previous studies based on L. monocytogenes strain ScottA and LO28 barotolerant isolates showed that there is a phenotypic and genotypic variation between barotolerant isolates [36,37,43], which indicates that there are variety of factors which can lead to pressure tolerance. Therefore, it is a challenge to link a genomic profile with pressure tolerance in different strains. Nevertheless, pan-genome association analysis results suggest that there was a genotypic difference between barotolerant and barosensitive strains. A total of 13 genes were significantly (p<0.01) associated with barotolerant strains. Interestingly, most of the barotolerant strain associated genes (Table S6) were located in the prophage regions in the studied genomes. In addition, extension search of these genes within all complete L. monocytogenes strain genomes showed that some phage genes (OCPFDLNE_02579 and OCPFDLNE_02580) were only seen in strains that contain inactivated CRISPR/Cas system (Table S6). It has been shown that prophages can provide increased biofilm formation and can be beneficial for stress coping in both L. monocytogenes and Escherichia coli [44,45]. We also observed that barotolerant strains harbored a slightly higher number of prophages compared to barosensitive strains (based on average number of prophages in barotolerant and barosensitive strains). Therefore, we believe that prophages might play a role in barotolerance in L. monocytogenes.
It was also interesting to see that two or more CRISPR-Cas systems, self-targeting spacers, and anti-CRISPR genes were seen together only in the three most barotolerant strains (RO15, MB5, and 2HF33). Since the existence of a self-targeting spacer indicates that the activity of the CRISPR/Cas systems is inhibited [46], it can be predicted that anti-CRISPR proteins detected in barotolerant strains inactivate CRISPR/Cas system. Previous studies have also shown that anti-CRISPR proteins are able to inhibit several types of CRISPR systems [33,47]. The observed active transcription of anti-CRISPR genes in RO15 prophages based on RNA-seq may also support this prediction. It has been shown that bacteria tend to lose CRISPR/Cas systems, since they can be harmful to the host itself. Immunopathological effects through partially matching spacers were reported unless the phage carries anti-CRISPR genes [48]. The anti-CRISPR genes found in prophages and inhibition of the CRISPR/Cas system may be the reason why barotolerant strains carry type I or type II CRISPR/Cas systems. We observed that more than 50% of the L. monocytogenes strains that harboring type I or type II CRISPR/Cas system also had anti-CRISPR genes (Table S6). CRISPR/Cas systems can have alternative roles in bacteria, such as roles in pathogenesis or role in modulating biofilm formation [49], however the role of inhibited CRISPR/Cas system is not known. More studies are required to support the relation between inhibited CRISPR/Cas systems and barotolerance.
RNA-seq revealed the general transcriptome profile of strains RO15 and ScottA. Significantly upregulated heatshock genes after pressure treatment (chaperone and clp protease genes) in ScottA (Table S2) might indicate that proper folding of proteins plays an important role in cell recovery after pressure treatment. A previous study also showed that the expression levels of Clp proteins were significantly higher in mutant HPP-tolerant ScottA isolate compared to wild type ScottA [43].
A previous study concluded that antibiotic resistant L. monocytogenes strains are more tolerant to pressure at 400 MPa [19]. Here, annotation results did not show any strain-specific antibiotic resistance gene. Nevertheless, multiple sequence alignment of predicted antibiotic resistant genes showed that there were slight differences in amino acid sequences across the strains ( Figure S3). However, it is not known that these amino acid differences cause antibiotic resistance advantage, and further studies are required to link the amino acid sequence variations of antibiotic resistance genes and pressure tolerance.

Conclusions
In this study we noticed that barotolerance is manifested at pressures lower than 600MPa. Strain RO15 was identified as the most barotolerant strain for 400 MPa 1 minute pressure treatment. Genome sequence of seven new strains and genome comparison of 10 strains revealed that the three most barotolerant strains have CRISPR-Cas genes and anti-CRISPR genes in their genomes. In addition, the average number of prophages was slightly higher in barotolerant strains compared to barosensitive strains. Furthermore, we have predicted 10 phage genes that might be related to pressure tolerance based on the pan-genome association test. Therefore, we conclude that prophages and inhibited prophage defense systems may be linked to pressure tolerance. The observation described above may argue for the use of phage cocktails as a pretreatment before HPP and will allow to reduce pressure and holding time.

Strains, Growth conditions and Pressure Treatment
Strains Scott A (CIP103575) and EGD-e (CIP107776) were obtained from Centre de Ressources Biologiques de l'Institut Pasteur, Paris, France, and strain F2365 (LMG23356) from Laboratorium voor Microbiologie, UGent, Gent, Belgium. Strains RO4, RO15, AB120 and AB199 are from Dunarea de Jos University of Galati, Romania, and have been isolated in the Promise FP7 project either from food illegally introduced to Romania or from meat processing environments. All other strains were isolated from the environment as indicated in Table 2. The strains were stored on Microbank beads with cryopreservatives (Pro-Lab Diagnostics, Toronto, Canada) at -80 °C prior to use. The bacteria were passaged twice in tryptic soy broth (TSB) supplemented with 0.6% (w/v) yeast extract (TSBYE; Oxoid, Basingstoke, Hampshire, England). Prior to pressure treatments, bacteria were grown overnight in 50 mL TSBYE at 37 °C with rotary agitation (150 rpm), resulting in cells in early stationary phase and with a target concentration of approx. 10 9 cfu/mL. Aliquots (10 mL) of the cultures were packaged in sousvide plastic pouches and sealed without using vacuum, and 10 mL was transferred to 15 mL falcon tubes to serve as untreated controls.
HPP was carried out using the QFP 2L-700 (Avure Technologies Inc., Columbus, USA). The cylindrical pressure vessel had 10 x 25.4 cm dimensions, 2 L capacity and 690 MPa upper pressure limit. Temperature of the pressure medium (water) was tracked with a K-type thermocouple located on the external surface of the samples. A holding time of 1 minute was used, pressures of 400 and 600 MPa, and ambient vessel water temperature (20-22 °C). Due to adiabatic heating, water temperatures in the middle of the vessel at the end of pressure treatment had risen to 31-33 °C after the 400 MPa treatment, and 36-38 °C after the 600 MPa treatment. Straight after pressure treatment, 400 MPa pressurized and untreated samples were serially ten-fold diluted in TSBYE and plated in triplicate on tryptic soy agar with 0.6% yeast extract (TSAYE; Oxoid) by using a spiral plater (Eddy Jet; IUL Instruments, Barcelona Spain). Pressurized samples at 600 MPa were additionally plated manually (100 μL) without being diluted. TSAYE plates were incubated at 37 °C for 48 hours prior to counting the colonies and estimating bacterial inactivation. Two consecutive trials of the methodology were performed (Exp. 1 and 2).

DNA extraction
DNA of each strain was extracted from 5 ml of a culture grown overnight in BHI broth at 37 °C with aeration on a rotary shaker. Wizard Genomic DNA purification kit (Promega, Madison, WI, USA) was used according to the manufacturer's instructions.
Library preparation, genome sequencing, de novo assembly, base modification detection and motif analysis L. monocytogenes strain RO15 was sequenced using PacBio RSII (Pacific Bioscience, Menlo Park, CA, USA) and Illumina Miseq (Illumina, San Diego, CA, USA). Genomic DNA PacBio library was prepared according to the manufacturer's Template Preparation and Sequencing Guide using DNA Template Prep Kit 2.0 and DNA/Polymerase Binding Kit P6. The Nextera XT protocol was used for the enzymatic genomic DNA fragmentation and Illumina sequencing library preparation. Pacbio reads were assembled using HGAP3 protocol [50] in SMRTPortal 2.3.0. Obtained assembly of chromosomal and phage sequences were checked and circularized using Gap4 program [51] and finally the chromosomal DNA sequence was set to start from dnaA gene. cutadapt v1.8.1 [52] was used with -m 200 and -q 25 options to trim Nextera adapter sequences and quality filtering. Trimmed reads were mapped against circularized chromosomal and phage sequences using bwa-mem [53]. Short indels within homopolymeric regions were corrected using pilon v1.16 [54]. Average sequencing coverages at the whole genome level were 308X in Pacbio data and 157X in Illumina data, respectively. DNA base modifications and motifs were analyzed using Modification and Motif analysis protocol implemented in SMRTPortal 2.3.0 with default parameters.

Pressure Treatment for RNA-seq experiment and RNA extraction
For both strain RO15 and ScottA, 10 ml of cells were grown in TSBYE until early stationary phase (each triplicated with new cultures). Samples were packed in sterile sous-vide pouches and packaged without applying vacuum. They were put in the refrigerator 30 minutes before pressure treatment. After treatment at 200 and 400 MPa for two, eight and sixty minutes at 20-22°C, samples were transferred to sterile 15 ml Falcon tubes. From each sample, 0.5 ml were transferred to a 2 ml Eppendorf tube prefilled with 1 ml RNA Protect (Qiagen, Hilden, Germany). They were vortexed, incubated for 5 minutes at room temperature, and frozen at -80 °C. Cells previously stored in RNA Protect (Qiagen) were pelleted by centrifugation for 10 minutes at 5000 g. RNA extraction was performed with NucleoSpin RNA kit (Macherey-Nagel, Düren, Germany) according to manufacturer's instructions with some modifications in the cell disruption phase. The cell pellets were suspended with 700 µl RA1 buffer and 7 µl β-mercaptoethanol (Sigma-Aldrich, Saint Louis, MO, USA). Cells were then mechanically disrupted using Lysing Matrix B tubes (MP Biomedicals) and FastPrep 24 tissue homogenizer (MP Biomedicals, Irvine, CA, USA) at 6 m/s for 3 x 30 sec. Cells were rested on ice for five minutes between cycles. After spinning the cells briefly and transferring the supernatant to NucleoSpin filter column, manufacturer's protocol was followed. Quantity and quality of RNA extractions were analyzed using the Agilent 2100 Bioanalyzer and RNA 6000 Nano kit (Agilent, Santa Clara, CA, USA).

RNA-seq library preparation
Ribosomal RNAs were removed from the total RNAs (ScottA 9.4 µl; RO15 14 µl) with Ribo-Zero rRNA Removal kit for bacteria (Illumina) according to manufacturer's instructions with 1/3 (ScottA) or 1/2 (RO15) volumes of kit solutions. rRNA-depleted RNA was purified with RNeasy MinElute Cleanup Kit (Qiagen) according to modified protocol of the Ribo-Zero kit manual and eluted in 12 µl of RNA-free water. Eight µl of rRNA-depleted RNA was used to prepare RNA sequencing libraries with SENSE Total RNA-seq Library Prep Kit for Illumina (Lexogen, Vienna, Austria) according to manufacturer's instructions. In reverse transcription and ligation phase, incubation time was extended to two hours. After second-strand synthesis, purification and size-selection of the libraries were performed with 13 µl of Bead Diluent and 27 µl of Purification Solution. PCR amplification program was slightly modified by increasing the denaturation times: from 30 to 60 sec in the beginning, and from 10 to 30 sec during the cycles; cycle number was increased to 40

Genome annotation, prophage prediction and genome alignment
Assembled genomes were annotated using Prokka v1.13 [63] with default options. To improve functional annotation PANNZER2 annotation web server [64] was used. CRISPR repeat regions were detected using CRISPRCasFinder version 4.2.19 [65]. Anti-CRISPR genes were annotated using BLAST alignment against known anti-CRISPR genes [33]. The prophage prediction for all genomes was done using PHASTER web tool [66]. Serotype of the selected strains were predicted by checking the marker primers [24] in the genome using EMBOSS primersearch v 6.6.0 [67]. Multilocus sequence typing (MLST) based on 7 loci (MLST), clonal complex (CC), core genome MLST (cgMLST) and lineages of strains were assigned based on MLST and the cgMLST schemes developed by Moura et al. [26,68] by uploading the sequences to BIGSdb-Lm webserver (https://bigsdb.pasteur.fr/listeria/listeria.html). Antibiotic resistance genes were detected using Resistance Gene Identifier (RGI) tool with CARD database [69]. Heuristic neighbor-joining phylogeny tree based on concatenated sequences of the MLST gene fragments was created using FastTree v2. 1 [70]. Whole genome alignment and ring figure was created using CGView Comparison Tool [71]. The ANIb scores were calculated using JSpeciesWS webtool [72].

Pan-genome and Pan-GWAS analysis
The pan-genome analysis was done using Roary Pan-genome Pipeline [73] with default settings, to create alignment of core genes using PRANK [74] the "-e" setting was used. The core genes alignment was used for phylogenetic tree construction using FastTree [70]. Pan-genome-wide association analysis was done using Scoary [75] with default settings.

Availability of data and material
All sequencing data and assembled genomes have been deposited in the European Nucleotide Archive (ENA) under accession code PRJEB35939. Strains RO4, RO15, AB120 and AB199 are from Dunarea de Jos University of Galati, Romania, and are available upon request to anca.nicolau@ugal.ro. Strains C7, MB5, and 2HF33 are from Nofima -Norwegian Institute of Food, and are available upon request to trond.lovdal@nofima.no.    Figure S3. Table S4 shows identified antibiotic resistance genes and percentage identity of matching region in all strains using CARD database. Figure S3 shows average distance trees from alignments of protein sequences of identified antibiotic resistance genes.
Additional File 6: Table S5. T-test against RO15 log 10 reduction. Additional File 7: Table S6. Barotolerant strains specific genes based on pan-genome wide association analysis. Table shows genes that have p < 0.01 based on pan-genome wide association analysis for barotolerant strains. Gene IDs of strain RO15 was used for this table. Orthologs of these genes with 95% identity were also seen in strains MB5 and 2HF33. In addition, the list of other publicly available L. monocytogenes genomes, where these genes were seen, are also provided. Figure 1 Viable cell counts bar chart. Viable cell counts log10(cfu/ml) of untreated controls and samples treated for 1 minute at 400 MPa (a) or 600 MPa (b). Gray bar represents control samples, blue bar represents treated sample. Error bars represent standard deviation (ANOVA; *, p < 0.02).

Figure 2
Circular map of the L. monocytogenes strain RO15 chromosome with predicted prophage regions, CRISPR region and anti-CRISPR genes. Five prophage regions were predicted in the chromosome of L. monocytogenes strain RO15 using PHASTER tool, shown in figure with blue and green color. Two anti-CRISPR gene regions were annotated within the prophage region 4 and 5. RliB-CRISPR/Cas system was located at position 525-527 kbp, following CRISPR-I system located the position 535-547 kbp with cas1-8 genes and a CRISPR array with 54 spacers.

Figure 3
Genome comparison of the selected nine strains mapped against strain RO15. The genomes were aligned using BLAST and visualized as ring figures. Each ring represents the genome alignment to strain RO15. From inside to outside, strains are located according to increasing hits of high percent identity to strain RO15 (Order from inside to outside: strain F2365, ScottA, AB199, AB120, RO4, MB5, EGD-e, C7, and 2HF33. Color represents BLAST hit identity percentage).

Figure 4
Visualization of CRISPR systems in strain RO15, 2HF33, and MB5. The figure shows a representation of CRISPR systems in strain RO15, 2HF33, and MB5. These three strains contained more than one CRISPR system, while the rest of the selected strains in the study contained only the RliB-CRISPR system. Arrows represent genes, rectangles represent repeat/spacer arrays, the numbers below the rectangles indicate number of spacers in the array. For simplification, sizes of the arrows do not correspond to the actual size of the genes.

Figure 5
Pan-genome comparison of selected strains. A phylogenetic tree was created based on core genome alignments of the selected strains. The matrix shows presence (blue) and absence (white) of core and accessory genes. Core genes are found in all strains, accessory genes are found in at least one strain.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.