HPP-induced transcriptome response during recovery of Listeria monocytogenes

High-pressure processing (HPP) is a commonly used technique in the food industry to inactivate pathogens, including L. monocytogenes. It has been shown that L. monocytogenes is able to recover from HPP injuries and can start to grow again during long-term cold storage. To date, the gene expression proling of L. monocytogenes during HPP damage recovery at cooling temperature has not been studied. In order identify key genes that play a role in recovery of the damage caused by HPP treatment, we performed RNA-sequencing (RNA-seq) for two L. monocytogenes strains (barotolerant RO15 and barosensitive ScottA) at nine selected time points (up to 48 h) after treatment with two pressure levels (200 and 400 MPa).


Introduction
Listeria monocytogenes is a foodborne bacterial pathogen that poses a particular challenge to the food industry due to its ubiquitous nature and capability of adapting to various inhospitable environmental conditions related to food matrices and food processing environments [1,2,3]. Transmission of this bacterium to humans generally occurs via consumption of contaminated food, especially ready-to-eat (RTE) foods that do not undergo thermal treatment during the manufacturing process, such as RTE salads, dairy products from raw milk, vegetables, and fruits. Another foodproduct category of that sliced and packed meat products may be contaminated with L. monocytogenes is [4]. L. monocytogenes can cause listeriosis, a disease associated with a high number of hospitalization cases and mortality rates of 20-30% among people with weakened immune systems [4]. L. monocytogenes can resist a wide range of environmental conditions [5] and its ability to grow at refrigeration temperatures increases the risk of listeriosis [6].
The increasing demand of consumers for minimally processed foods, with fresh-like sensorial and nutritional properties, requires the implementation of alternative food processing techniques such as high-pressure processing (HPP). HPP is a relatively new, non-thermal processing technique that shows remarkable results with respect to pathogen inactivation and minimum impact on food quality [7,8].
It has been reported that HPP causes morphological, structural, physiological, and genetic changes or damages to L. monocytogenes cells [9]. However, the susceptibility of L. monocytogenes to HPP varies between growth phase [10], strength of the cellular envelope [11], genomic features [12], and individual strains [13]. In addition, food matrix type, temperature, water activity, compression and decompression rates, applied pressure and holding time, and other extrinsic factors have an impact on inactivation of L. monocytogenes cells by HPP [9]. Several studies reported the potential of sublethally injured L. monocytogenes cells to recover after HPP and grow within the storage period even under refrigeration conditions [14,15,16,17,18]. Bozoglu et al. [19] showed that even though sublethally injured bacteria could not be detected immediately after HPP treatment of up to 550 MPa. However, injured but viable cells may be present in the pressurised samples as the authors detected growth after 6 days at 4 °C and already after 1 day at 22 °C and 30 °C.
Therefore, for an e cient decontamination process, the additional hurdles to increase e ciency of HPP and/or to prevent outgrowth of sublethally injured bacteria should be considered. In this context, it may be of interest to treat L. monocytogenes cells with antimicrobial agents that compromise cell wall and/or membrane and thereby render bacteria more sensitive to HPP and inhibit recovery. Such antimicrobial agents may include bacteriocins [20], essential oils [21,22,23], plant extracts [24], bacteriophages [25,26], lysozyme [27,28], lactoferrin [29], and lactoperoxidase [30].
Gene expression pro ling of the response of L. monocytogenes to HPP has previously been studied by RNA-seq [12], microarray [31,32], and qPCR [33]. For example, it has been shown that a mutation in ctsR causes barotolerance and a ctsR deletion mutant of L. monocytogenes shows increased expression of Clp protease and PTS system genes after HPP [31]. Similarly, we previously reported that heat-shock and Clp protease family genes were upregulated after HPP [12]. In contrast, Bowman et al. [32] reported downregulation of heat-shock and PTS system genes after HPP. The previous studies used relatively higher temperatures for storage after HPP (⩾ 15 °C) compared to common coldstorage applications in the food industry. We recently showed genetic differences between barotolerant and barosensitive L. monocytogenes strains, which may explain their different HPP sensitivity [12]. Hence, in the present study, we investigate the transcriptional response to HPP and the differences in gene expression pro les between barotolerant and barosensitive L. monocytogenes strains during recovery. We selected L. monocytogenes strains RO15 (barotolerant) and ScottA (barosensitive) that were already analysed previously [12]. This is the rst study to perform time-series RNA-seq using both barosensitive and barotolerant strains monitoring gene expression pro les during recovery of an HPP insult at 8 °C. We aimed to identify candidate genes that would be involved in the recovery of L. monocytogenes after HPP treatment.

Log reduction testing of the strain RO15 and ScottA
We previously reported that strain RO15 is more resistant to treatment of 400 MPa for 1 minute compared to several other L. monocytogenes strains including strain ScottA [12]. Here, we sought to test the susceptibility of both strains to pressure treatment at 200 and 400 MPa for 8 minutes at 8 °C. While a treatment with 200 MPa was ineffective for inactivation of both strains, 400 MPa signi cantly reduced log 10 (cfu/ml) for both RO15 (5.78) and ScottA (7.04) compared to untreated samples (p < 0.05). The log reduction difference between the two strains was also statistically signi cant (p < 0.05; Fig. 1).

Differential Expression Analysis
After the HPP treatments, samples were taken at nine time points (0, 5, 10, 30, 45, 60 min and 6, 24, 48 h) and RNAseq performed. Principal component analysis (PCA) of per gene read count data showed that there was a clear separation between HPP-treated samples and control samples for both 200 MPa and 400 MPa ( Figure S1). In addition, we also saw clustering between early time points (0, 5, and 10 min), middle time points (30,45, and 60 min), and late time points (6,24, and 48 h) for 200 MPa treatment on a PCA plot ( Figure S1).
Pairwise differential expression analysis between the treated samples and corresponding controls at all time points showed that a large number of genes were signi cantly differentially expressed (p-value ≤ 0.05) after HPP ( Figure S2, Table S1, S2, S3, S4).

Differentially expressed genes and GO enrichment analysis
To gain a general perception of the functional groups of the differentially expressed genes, GO enrichment was performed (Fig. 2, Table S5). The most signi cantly enriched GO terms for upregulated genes were cobalamin biosynthetic process, divalent inorganic cation transport, and transition metal ion transport for both strains (Fig. 2). These GO terms were enriched at several time points in both strains after 200 and 400 MPa treatment. The main upregulated genes responsible for these GO terms enrichment were found in a large operon (OCPFDLNE_01234 -OCPFDLNE_01251 in the RO15 strain, LMOSA_20560 -LMOSA_20730 in the ScottA strain), including cobalamin biosynthesis genes, Cobalt ABC transporter, and Cobalt transport-related genes ( Figure S3).
In L. monocytogenes RO15 HPP-upregulated genes were enriched at most time points in GO terms "phosphoenolpyruvate-dependent sugar phosphotransferase system (PTS system)" and "carbohydrate transmembrane transport" both after the 200 and 400 MPa treatment (Fig. 2, Table S5). In detail, upregulation was observed for the genes for fructose-, mannose-, galactitol-, cellobiose-, and ascorbate-speci c PTS systems. Enrichment of these GO terms was also seen in the HPP-treated samples of L. monocytogenes ScottA strain taken, however, at later time points after HPP (both 200 and 400 MPa) (Fig. 2).
HPP caused upregulation of protein folding, chaperone, and peptidases genes, such as clpE, clpP, groEL, groES, hrcA, dnaK, and dnaJ, at 200 MPa at almost all time points and at 400 MPa at the early time points as re ected by an enrichment in the GO term "protein folding" (Fig. 2, 3c, 4c).
Signi cant downregulation was observed for cell division genes divIC, dicIVA, ftsE, and ftsX ( Fig. 3i, 4i). In addition, in downregulated genes we observed a signi cant (p < 0.05) enrichment of the GO term "ATPase activity" (Table S5) at almost all timepoints for both pressure levels.
In addition to enriched GO terms, RNA-seq data was also analysed for speci c gene families such as transcription factors ( Fig. 3b and 4b, Figure S4) and transcription machinery genes ( Fig. 3f and 4f, Figure S5). The gene hrcA (encoding heat-inducible transcription repressor HrcA) was upregulated at all time points for both RO15 and ScottA strains after 200 MPa, while after 400 MPa treatment upregulation was seen only in ScottA for the early time points (Fig. 3b, 4b). We also observed that gene prfA encoding the master regulator of virulence genes in L. monocytogenes was upregulated in both strains after both pressure treatments (Fig. 3b, 4b). Interestingly, only one of the manR genes encoding the transcriptional regulator ManR, which is found in two different copies in the genomes of ScottA and RO15, was upregulated in RO15 strain after HPP treatment but not in ScottA ( Fig. 3b and 4b). With respect to the transcription machinery genes, rpoD encoding the RNA polymerase sigma factor RpoD was upregulated in both strains after HPP ( Fig. 3f and 4f).
We also focused on mechanosensitive channel genes. The gene encoding a putative mechanosensitive channel gene of large conductance (mscL) was upregulated at early time points after 200 MPa treatment in RO15, while this upregulation was not seen in ScottA (Fig. 3h). After 400 MPa, both strains had similar mscL expression patterns with upregulation at late time points (Fig. 4h). The homologue for a mechanosensitive channel of small conductance (ykuT) was only upregulated in RO15 at 48 h after 200 MPa treatment (Fig. 3h).
To see the difference between the responses to the different pressure levels, we identi ed genes that were upregulated after 200 MPa treatment but not after 400 MPa and vice versa for each time point. In both strains, the genes that were upregulated after 200 MPa but not after 400 MPa at early time points were mainly related to translation (Table S6, Table S7). Interestingly, translation-related genes were upregulated after 400 MPa but not after 200 MPa in the RO15 on late time points (Table S6, Table S7). We observed upregulation of hpf gene (encoding ribosome hibernation promoting factor) in ScottA even at the time point 48 h (Fig. 4j). In addition, we also observed several cobalamin biosynthesis and PTS-related genes were upregulated at early time points after 400 MPa but not after 200 MPa in both strains (Table S6, Table S7).
Genes without orthologs within both strains were mainly phage genes and hypothetical genes. In both strains, phage genes were mostly upregulated after HPP ( Figure S6, Figure S7). We previously reported that barotolerant strains harbour both CRISPR-Cas systems and anti-CRISPR genes [12]. However, upregulation of Cas genes was observed in neither of the two strains whereas anti-CRISPR genes (acrIIA1 and acrIIA2) were signi cantly upregulated after HPP in RO15 ( Figure S6).
Non-coding RNA (ncRNA) RNA-seq read coverage plots showed that a very large amount of RNA-seq reads were mapped to non-coding regions, especially for RO15. Further examination showed that, on average, ~ 53% (ranging from 21-86%) of all RNA-seq reads for RO15 samples were mapped to the small non-coding RNA (ncRNA) Rli47. Similarly, ~ 28% (ranging from 6-72%) of the RNA-seq reads in samples of ScottA mapped to Rli47 (Table S8, Table S9). Thus, we additionally performed expression analysis of ncRNAs. We observed that Rli47 transcript levels were upregulated in response to pressure treatment in both strains ( Figure S8). Similarly, levels of LhrA ncRNA were upregulated in both strains at the early timepoints. Interestingly, expression of Rli53 was upregulated in RO15 after the pressure treatment, while no upregulation was seen in ScottA.
Gene regulatory networks based on RNA-seq data One of our goals was to understand the regulatory networks involved in the response to HPP treatment in L. monocytogenes strains, RO15 and ScottA. Consensus gene network was created using the time-series expression data for all differentially expressed genes in both strains. This resulted in a total of 3661 gene network links (1506 genes and 3661 edges) for strain RO15 and 3427 gene network links (1389 genes and 3427 edges) for strain ScottA (Table S10). Interactive visualizations can be seen on https://icemduru.github.io/RO15_gene_network and https://icemduru.github.io/ScottA_gene_network. Moreover, we clustered the genes based on the network data (Table S10) using network clustering algorithm Map equation [34]. For RO15, 151 clusters were predicted in the gene network (Table S11), while for ScottA, 128 clusters were predicted (Table S12).
For both strains, heat shock and chaperone-related genes were clustered together (Cl6 in RO15 and Cl9 in ScottA, Fig. 5, 6). De novo motif discovery analysis resulted in a number of signi cant motifs (E-value < 0.05) for the upstream regions of heat shock clusters (Cl6 in RO15 and Cl9 in ScottA), and one of the motif was signi cantly (Evalue < 0.05) similar to the CtsR motif (Table S13, Table S14) from the PRODORIC database [35] in both strains. This indicates CtsR is a regulator for protein-folding genes in these strains. We also observed that ctsR was linked to heatshock genes based on gene network inference (Table S11, Table S12). Notably, nagA and nagB are placed in the heat shock cluster (Cl9) in ScottA providing a hint that co-expression of protein folding genes and peptidoglycan biosynthesis genes after the pressure treatment was required together for recovery in ScottA. In addition, we observed that the heat-shock cluster (Cl9) in ScottA was linked to Cl4 (Fig. 6), which includes stress-related genes.
The prophage genes were highly interconnected in both strains; almost all genes of the three phages found in the same cluster in ScottA (Cl2, Fig. 6) indicate all three prophages show similar gene expression reactions in ScottA. Similarly, prophage 3 and prophage 5 genes were highly linked in RO15 (Cl2 and Cl3, Fig. 5). The prophage cluster (Cl2) was also linked to Cl4 (Fig. 6) in ScottA, which contains universal stress protein UspA genes (uspA), indicating that phage induction was connected to stress response in ScottA. The prophage cluster (Cl3) was linked to Cl5 in RO15 (Fig. 5), which includes mscL, i.e. the gene for a mechanosensitive channel gene of large conductance, and cspA encoding a cold-shock protein.
The genes of Cl9 in RO15 were enriched for The GO term "response to antibiotic". This cluster also included genes uspA (for universal stress protein A), virulence factor prfA (the master regulator of virulence in L. monocytogenes), and hpf (ribosome hibernation promoting factor), which all were signi cantly upregulated in both strains after HPP treatment. A similar cluster containing prfA and hpf was seen in ScottA (Cl6). This implies that stress response, virulence and ribosome hibernation are linked to each other and co-occurred with HPP treatment in both strains. For ScottA, de novo motif discovery for the upstream region of Cl6 genes resulted in two signi cant motifs (E-value < 0.05; Table S14), one of them being signi cantly (E-value < 0.05) similar to SigB motif from PRODORIC database [35]. ddPCR for validation for RNA-seq Results ddPCR was performed for the same ScottA samples that were used for RNA-seq to validate the differential expression results obtained using RNA-seq data. Nine genes were selected (Table S15)  The in uence of HPP on mutant strains` viability Based on transcriptional analysis, we decided to perform additional experiments on HPP resistance with the mutants carrying deletions in selected candidate genes. These candidate genes were selected based on the gene-expression data and a potential role in barotolerance and recovery from HPP damage. The selected candidate genes were: ykuT, encoding a putative mechanosensitive ion channel of small conductance (OCPFDLNE_01076 and LMOSA_19040 in RO15 and ScottA); and pbp2A, encoding for the penicillin-binding protein A2 (OCPFDLNE_02389 and LMOSA_2530 in RO15 and ScottA). For both genes, deletion mutants were generated in the widely used L. monocytogenes model strain EGDe (lmo1013 and lmo2229 in EGDe, respectively). Successful deletion was validated by genome sequencing ( Figure S11). In order to test the resistance of Listeria monocytogenes EGDe mutants to HPP, they were grown to stationary growth phase and subjected to HPP at 300 and 350 MPa with a holding time of 5 min. Susceptibility of the mutants to these treatments was assessed by detecting colony forming units. This revealed a reduction of ~ 5 log CFU/mL for Δlmo2229 cells and ~ 3 log CFU/mL for both Δlmo1013 and wild-type EGDe by the 350 MPa treatment and reduction of ~ 1 log CFU/mL for all three by 300 MPa treatment (Fig. 8). For both 300 MPa and 350 MPa treatments, we observed a signi cant (p < 0.05) difference in colony forming units between wild-type and Δlmo2229. No signi cant difference was observed between wild-type and Δlmo1013.

Discussion
High pressure processing (HPP) is commonly used in the food industry to inactivate foodborne pathogens and spoilage microorganisms. However, it has been reported that L. monocytogenes is able to recover after HPP treatment during long-term storage [18,19,36]. To study gene expression response during early recovery of L. monocytogenes and identify genes that are important for recovery, we performed RNA-seq of samples taken at different time points after HPP at two different pressure levels (200 MPa and 400 MPa). To account for strain-dependent variation in the HPP response, experiments were performed with two strains: RO15, a strain that was shown to be more resistant to HPP than others; and ScottA, a strain that is more sensitive to HPP [12].
L. monocytogenes can recover within one day after injury caused by HPP for 10 min at 450 MPa and 45 °C [19]. Recovery of L. monocytogenes was also observed after six days of storage at 4 °C for injured bacteria following treatment with HPP at 550 MPa at 45 °C for 10 min [19]. Our results show that viability was unaffected in the two strains tested after HPP treatment at 200 MPa and 8 °C for 8 ( Fig. 1), but a marked reduction in viable counts of several log 10 units was observed for both strains treated at 400 MPa. However, the reduction in viable counts was signi cantly higher in ScottA compared to RO15, which supports our previous observation that RO15 is more pressure tolerant [12].
The time-resolved RNA-seq data allowed us to perform gene-network analysis. To summarize the gene networks, we clustered genes assuming that genes within a cluster, and in linked clusters, are functionally related or interact during recovery from HPP [37,38,39]. HPP mainly affects expression of protein folding genes, PTS system genes, prophage genes, and cobalamin biosynthesis genes. We saw that in both strains, stress response genes, virulence genes, and ribosome hibernation promoting factor hpf gene were strongly linked to each other, indicating that during recovery from HPP, co-expression of these three factors was needed. It has also been reported that the general stress sigma factor B (σ B ) regulates hpf, prfA (encoding listeriolysin regulatory protein), and UspA1,2 (encoding universal stress proteins A1 and A2) [40]. Environmental stress activates σ B , which regulates more than 200 genes [40,41]. In line with this, based on de novo motif discovery analysis, SigB transcription factor binding site-like motif was found in the upstream regions of the gene cluster (Cl6; Fig. 6, Table S14), which includes prfA and hpf in ScottA. This indicates that a general stress response was activated in L. monocytogenes by σ B after HPP.
We have reported that a large portion (up to ~ 85%) of RNA-sequencing reads were mapped to Rli47 ncRNA, which was upregulated in both strains after pressure treatment. Similarly, previous studies have also reported that up to ~ 90% of all RNA-seq reads map to Rli47 ncRNA in L. monocytogenes [42,43]. It has been shown that Rli47 plays a role in the response to acid stress [43] and oxidative stress [44]. In line with these observations, our data suggests that Rli47 is also involved in HPP recovery based on high expression level after HPP treatment. It has also been shown that Rli47 is regulated by σ B [45]. This supports our observation of general stress response activation by σ B after HPP.
In addition, Rli53 expression was upregulated in RO15 but not in ScottA. Rli53 has been associated with antibiotic resistance [46]. Our results indicate that Rli53 may also play a role in pressure resistance in RO15.
Cobalamin biosynthesis was the most signi cantly enriched GO term for upregulated genes in both strains. It has been shown that cobalamin plays a protective role against oxidative stress in bacteria [47]. Cobalamin was also shown to be an essential cofactor for propanediol and ethanolamine utilization [48]. Signi cant downregulation of cobalamin biosynthesis genes in L. monocytogenes has been reported in response to Rli47 deletion [45]. Hence, upregulated cobalamin synthesis genes after HPP in this study can be related to increased levels of Rli47, which is regulated by σ B . We therefore predicted that cobalamin biosynthesis genes were upregulated as part of the general stress response of HPP.
Stress conditions have been shown to induce prophages in L. monocytogenes [49,50]. Upregulation of prophage genes in both strains after HPP indicates that pressure stress also induces prophages. In addition, co-regulation of different prophages within the same host has also been shown in L. monocytogenes [50]. Similarly, our gene-network inference suggests that prophages were linked to each other within strains, coexpression of prophages genes were observed. Based on pan-genome analysis, we previously proposed that prophages and anti-CRISPR genes may play a role in pressure resistance in L. monocytogenes [12]. In this study we observed that both anti-CRISPR genes (acrIIA1 (OCPFDLNE_02770, OCPFDLNE_02583) and acrIIA2 (OCPFDLNE_02582)) in RO15 were upregulated after HPP. GO enrichment analysis of upregulated genes indicated that PTS systems were activated in both strains for most of the time points during the recovery phase. Upregulation of PTS genes has also been reported for other stress conditions in L. monocytogenes based on transcriptome analysis [51,52]. Upregulated PTS systems were mostly galactitol-, fructose-, and mannose-speci c PTS systems. These carbon sources play a role in cell-wall biosynthesis [53,54]. Thus, upregulation of these sugar transporters may be an indication of increased uptake of these sugars for cell-wall biosynthesis and as a carbon source to allow recovery from injury caused by HPP.
CtsR is a negative regulator of heat-shock genes, mainly of the clp family of genes, and has been shown to be directly involved in barotolerance of L. monocytogenes [55,56,57,58,59]. Deletion of ctsR led to an increase in barotolerance of ScottA by 5 orders of magnitude [58]. In addition, upregulation of PTS, heat-shock, and clp family genes has been reported for a ScottA ctsR mutant [31]. Furthermore, ctsR is reportedly regulated by σ B in Bacillus subtilis [41]. Our results show upregulation of ctsR expression in both strains at some time points after treatment with 200 MPa but not with 400 MPa. Moreover, upregulation of genes was observed for heat-shock proteins of the clp family and chaperones in samples treated at both 200 MPa and 400 MPa in both strains (Fig. 3c, 4c). Especially clpE was one of the most signi cantly upregulated genes at several time points. It has been shown that heat-shock proteins are needed to deal with misfolded proteins, prevent cellular damage, and help cell recovery during pressure treatment [60]. Our observation that genes for heat-shock proteins are upregulated indicates a similar role in protection and recovery of L. monocytogenes to/after HPP treatment.
It has been shown that antibiotic resistant L. monocytogenes strains are more resistant to 400 MPa pressure treatment compared to antibiotic-susceptible strains [13]. Our previous pan-genome study [12] also showed that barotolerant strains have slightly different amino acid sequences for norB encoding a protein involved in resistance against quinolones. Interestingly, different strains showed variations in their expression of norB. Signi cant upregulation of norB was observed in barotolerant RO15 at several time points after 400 MPa treatment, including early time points. However, norB was only upregulated at the 24 h time point after 400 MPa treatment in barosensitive ScottA. This supports the observation that differences in antibiotic resistance genes might provide a different barotolerance level within L. monocytogenes strains.
Ribosome damage can lead to cell death after HPP [61]. Ribosome hibernation, that is dimerisation of 70S ribosomes leading to translationally inactive 100S particles, has been reported to occur as L. monocytogenes adapts to different stress conditions [62]. Ribosome hibernation involves the gene product of hpf (ribosome hibernation promoting factor) and upregulation of hpf was seen in both strains. Downregulation of hpf was observed at the 48 h time point after 200 MPa treatment only in RO15. In addition, in RO15, the GO term "translation" was also mainly enriched for upregulated genes at the late time points in RO15 at both 200 and 400 MPa (Fig. 2). However, no enrichment of GO term "translation" was observed for upregulated genes in ScottA at 400 MPa (Fig. 2). Collectively, this indicates that L. monocytogenes keeps the translation inactive by inducing Hpf-mediated ribosome hibernation for a certain time after HPP. Moreover, there are differences between the strains in how long this hibernation lasts. The barotolerant strain RO15 seems to reactivate translation faster than the sensitive strain ScottA.
Based on morphological and physiological characterization, the cellular wall or membrane are targets to improve e cacy of HPP to inactivate L. monocytogenes [11]. We observed that peptidoglycan-synthesis genes such as, murG, murC, murD, and pbp2A were upregulated in both strains after HPP. Upregulation of peptidoglycan-synthesis genes with simultaneous downregulation of cell-division genes indicates that an active cell-wall repair occurs in both strains after HPP. pbp2A encodes a penicillin-binding protein that was shown to contribute to β-lactam resistance and cell morphology in L. monocytogenes [63]. To further investigate the role of this gene in response to HPP, a mutant carrying a deletion of the corresponding gene (lmo2229) was generated in L. monocytogenes strain EGDe and tested for resistance against HPP. The results show that the lmo2229 mutant was signi cantly more sensitive to HPP than the parental wildtype strain (Fig. 8). This shows that pbp2A is an important factor in L. monocytogenes for recovery after HPP.
Peptidoglycan of bacteria consists of a backbone of alternating N-acetylglucosamine (GlcNAc) and N-acetylmuramic acid (MurNAc) units interconnected with peptide side-chains [64]. It has been shown that bacteria are able to recycle N-Acetylglucosamine from peptidoglycan using the proteins encoded by nagA and nagB [53,65]. The genomes of L. monocytogenes RO15 and ScottA each contain two copies of nagA and nagB, respectively, and one of the copies of each nagA and nagB are organized in an operon. Expression of these nagAB copies were upregulated in both strains after HPP at 200 and 400 MPa. The second copy of nagA was not differentially expressed in either strains. Interestingly, the second copy of nagB gene (OCPFDLNE_02454) was only upregulated in RO15. This difference might partly explain the barotolerance difference between strains. Increasing protein levels of the NagB are associated with increased growth rate in E. coli [66]. Thus, more e cient biosynthesis of cell-wall peptidoglycans due to higher NagB levels may contribute to the higher barotolerance of RO15.
HPP creates a mechanical force that may result in deformation of the membrane. Mechanosensitive channels were shown to respond to membrane stress and help bacteria to cope with this stress [67]. We were intrigued by the observation that the mscL gene encoding a MS channel protein of large conductance was upregulated after 400 MPa pressure treatment in both strains. In addition, ykuT (encoding small MS channel protein) was upregulated at the 200 MPa 48 h time-point in RO15. However, the obtained lmo1013 mutant showed no signi cant difference in susceptibility/resistance to HPP indicating that the small MS channel protein was not directly involved in pressure resistance or only has a minor effect.
Gene expression pro ling under pressure treatment in Listeria was studied previously by Bowman et al. [32] using L. monocytogenes strain S2542. Notably, a signi cant negative correlation was observed for log 2 FC results between the study by Bowman et al. [32] and our study. Our RNA-seq results were validated using ddPCR (Fig. 7) and they are consistent for two different strains under different pressure levels and several time points. We speculate that discordance between the results could potentially arise from different growth conditions or different methods for measuring gene-expression levels.
Overall, these ndings may lead to new approaches to improve HPP e cacy. For example, we observed that the mannose phosphotransferase system (Man-PTS) was upregulated after HPP treatment. Man-PTS is the receptor for class IIa bacteriocins, such as pediocin or garvicin [68,69,70]. Thus, increased expression of these receptors may provide an opportunity to pre-treat food with IIa bacteriocins, which may increase susceptibility to HPP. However, dltD upregulation in RO15 may lead to incorporation of more alanine residues [71], which increases the positive charge and consequently reduces a nity to cationic antimicrobials and bacteriocins. Interestingly, dltD was downregulated in ScottA indicating pre-treatment might be more effective for barosensitive strains. In addition, among peptidoglycan biosynthesis genes, deletion of pbp2A causes signi cant susceptibility to HPP. Hence new approaches could be sought by using peptidoglycan cross-linking.

Conclusion
Recovery and outgrowth of L. monocytogenes in food after HPP treatment is a serious problem for the food industry. The mechanism of recovery of L. monocytogenes after HPP has not been studied by genome-wide transcriptional pro ling. Understanding how bacteria recover from HPP injury may help the food industry to develop new strategies for better inactivation of food pathogens. Here we reported a very detailed gene expression response of L. monocytogenes during recovery from HPP treatment using two strains (barotolerant and barosensitive), several time points, and two different pressure levels. Protein folding, PTS system, universal stress response, and cobalamin biosynthesis were the main activated functions in response to HPP treatment in L. monocytogenes. We showed that ncRNAs may also play a role in HPP injury recovery. Based on our results, several genes involved in barotolerance and recovery from HPP injury were predicted. Deletion of pbp2A suggests that it plays a role in barotolerance of L. monocytogenes. Further reverse-genetics experiments are required to validate our predictions based on RNA-seq.

Methods
High pressure processing for inactivation of strain RO15 and ScottA Pressure treatment testing the log reduction was carried out using the QFP 2 L-700 (Avure Technologies Inc., Columbus, USA) as previously described [12], except that a holding time of 8 min, vessel water temp of 8 °C, and pressures of 200 and 400 MPa were used for strain RO15 (isolated as described previously [12]) and ScottA (CIP103575, obtained from Centre de Ressources Biologiques de l'Institut Pasteur, Paris, France). Immediately after pressure treatment, pressurized and untreated samples were serially ten-fold diluted in tryptic soy broth with 0.6% yeast extract (TSBYE; Oxoid, Basingstoke, Hampshire, England) and plated in triplicate on tryptic soy agar with 0.6% yeast extract (TSAYE; Oxoid, Basingstoke, Hampshire, England) by using a spiral plater (Eddy Jet; IUL Instruments, Barcelona, Spain). Pressurized samples at 400 MPa were additionally plated manually (100 uL) without being diluted. TSAYE plates were incubated at 37 °C for 48 h prior to counting the colonies and estimating bacterial inactivation.
High pressure processing for RNA-seq Volumes of 350 mL BHI broth (Oxoid, Basingstoke Hampshire, England) were inoculated with L. monocytogenes RO15, and Scott A (CIP103575, obtained from Centre de Ressources Biologiques de l'Institut Pasteur, Paris, France) strain, respectively, previously sub-cultured in the same growth medium at 37 °C, at an OD 600 value of ~ 0.  Figure S12). At each time point, in order to stabilize the RNA, both treated samples (5 replicates) and corresponding controls (4 replicates) were transferred in 4 mL of RNA protect reagent (Qiagen, Hilden, Germany), incubated at room temperature for 5 min, pelleted by centrifugation at 5000 rpm and stored at -80 °C, until RNA was extracted.
RNA-sequencing RNA-sequencing was performed for 216 samples from the barotolerance experiment including three replicate samples for each treatment and time point (Table S16). Prior to RNA-seq library preparation, rRNAs were removed by hybridizing extracted RNA with DNA oligos complementary to 16S, 23S, and 5S rRNAs followed by digestion of resulting DNA-RNA hybrid molecules. Hybridization reactions included 1 µg RNA, 1 x buffer (16.7 mM Tris-HCl, pH8, 33.3 mM KCl), and 250 nM pooled DNA oligos in a total volume of 12 µl. Reactions were incubated at 95 °C for 2 min, slowly (0.1 °C/sec) cooled to 22 °C, incubated at 22 °C for a further 5 min, and placed on ice. For digestion of rRNAs, 1.5 µl of 10 x RNaseH buffer, 0.2 µl RNAseH (Thermo Scienti c, Waltham, Massachusetts, United States), 0.5 µl Ribolock RNase inhibitor (Thermo Scienti c), and 0.8 µl of water were added. Digestion reactions were incubated at 37 °C for 60 min, and inactivated at 65 °C for 20 min. Remaining unhybridized DNA oligos were removed with RapidOut DNA Removal kit (Thermo Scienti c) according to manufacturer's instructions. RNA-seq libraries were prepared using QIAseq stranded Total RNA Lib kit (Qiagen) using reaction volumes. Due to the large number of samples, puri cation steps were performed with a Magnatrix1200 pipetting robot (Magnetic Biosolutions) using precipitation on Dynabeads MyOne Carboxylic acid beads (Invitrogen, Carlsbad, California, United States), and 10% PEG after reverse transcriptase and second strand synthesis steps, and 9.5% PEG after the strand-speci c ligation step. Instead of the kit's adapter, a truncated TruSeq adapter was used in the strand-speci c ligation. Libraries were ampli ed using half of the puri ed ligation product (10 µl) in 1 x HF buffer with 0.2 mM dNTPs, 0.6 µM dual-index primer [72], and one unit of Phusion Hot Start II High-Fidelity DNA polymerase (Thermo Fisher Scienti c) in a total volume of 50 µl. The PCR protocol that was used was 98 °C, 30 s; 20 × (98 °C, 10 s; 65 °C, 30 s; 72 °C, 10 s); 72 °C, 5 min. Concentrations of ampli ed libraries were measured with a Qubit uorometer and dsDNA HS assay kit (Invitrogen, Carlsbad, California, United States), and size distributions visualized with Fragment Analyzer and High Sensitivity NGS Fragment Analysis kit (Advanced Analytical, Parkersburg, WV, USA). Ampli ed libraries were pooled into two batches, each including 108 samples. The rst pool was concentrated using an Amicon Ultra 100K column (Millipore, Burlington, MA, USA), puri ed once with 0.9 x AMPure XP beads (Beckman Coulter, Brea, CA, USA), and once with PEG (8-8.5%)/NaCl precipitation on Dynabeads MyOne™ carboxylic beads (Invitrogen). The second pool was puri ed twice with 0.9 x AMPure XP beads (Beckman Coulter). For both pools, size selection of 300-600 bp fragments was performed using BluePippin and 2% agarose gel cassette (Sage Science). Finally, the pools were puri ed with S-400 Microspin HR columns (GE Healthcare, Chicago, Illinois, United States). NextSeq 500 (Illumina, San Diego, CA, USA) was used to sequence the RNA-seq libraries twice. Altogether, seven sequencing runs were performed to produce 76 bp single-end reads.
RNA-seq data pre-processing and differential expression analysis.
RNA-seq reads obtained from all samples were processed using Trimmomatic v0.36 [73] to trim off the low-quality bases and lter out the adapter sequences. SortmeRNA v2.1b [74] was used to lter out ribosomal RNA reads. Reads were then mapped against the corresponding genomes (ScottA: GCA_000212455.1, RO15: GCA_902827145.1) using bowtie2 v2.3.4.3 [75] with default settings. Counting of reads per gene was performed using HTSeq v0.9.1 [76] with union mode. Hierarchical clustering of samples (HCA) based on Euclidean distances and principal-component analysis (PCA) of the samples was done for each L. monocytogenes strain as described in the manual for the DESeq2 R package v1.22.2 [77] on "regularized log" (rlog)-transformed read count data to visually explore sample relationships. One sample in ScottA (R_046) and two samples in RO15 (R_045, R_055) belonged neither to control nor to HPP clusters and were thus discarded as potential outliers. In addition, two RO15 samples (R_057, R_235) contained a very low number of CDS-mapped reads (approx. 0.05 million) and were also discarded (Table S17).
Differential expression analysis for each L. monocytogenes strain was performed using R package DESeq2 v1.22.2 [77]. The comparisons were made in multiple ways: 1) Comparing the HPP-treated samples to the corresponding controls for each time point and at different HPP level, separately; 2) Comparing the time dynamics of gene expression (as described in DESeq2 manual, http://master.bioconductor.org/packages/release/work ows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#timecourse-experiments) between control and treated samples, separately for 200 and 400 MPa series, to nd the genes that changed expression at least at one time point; the gene expression at time point 0 in control was taken as a proxy for gene expression before the HPP treatment (designated as T00). Genes were considered to be signi cantly differentially expressed if their adjusted p-value ≤ 0.05 and their log 2 fold change (log 2 FC) ≥ 0.6, (therefore, fold change FC ≥ 1.5). For the time dynamic comparison, genes were considered to be signi cantly differentially expressed at least at one time point if the adjusted p-value ≤ 0.05.

Gene regulatory network construction, clustering and visualization
The initial set of genes used to build the network contained 1,964 genes for strain RO15 and 1,852 genes for ScottA that changed expression between the corresponding controls and treatments at least at one time point (adjusted pvalue < 0.05). The gene expression values used to build a network were regularized log (rlog) transformed as described in DESeq2 [77] to stabilize the variance and normalize the count data. Eleven networks have been built using 11 publicly available algorithms (clr, genie3, aracne, pearson_corr, narromi, pcor, plsnet, tigress, llr-ensemble, elensemble) embedded into Seidr toolkit [78] and integrated into one network using Seidr. The hard threshold for the edge score was manually chosen to be 0.4 based on the criteria described here (https://seidr.readthedocs.io/en/latest/source/threshold/threshold.html) leaving only 1,506 genes, 3,661 edges for strain RO15, and 1,389 genes, 3,427 edges for strain ScottA in the nal weighted undirected network. Infomap v0.19.26 [34] two-level clustering with options "--bftree − 2 --ow-network -N 10000" was used to nd clusters (modules) of genes in the network and the ow between the modules (the strength of interactions between the modules). The two-level representation was visualized using the map&alluvial generator (http://www.mapequation.org/apps/MapGenerator.html). The network was visualized and centrality metrics calculated for the nodes (e.g. Degree, Betweenness and Closeness) using Cytoscape [79]. MEME suite v5.0.5 [80] was used for motif analysis. The upstream regions of genes ranging from 50 to 300 nucleotides were extracted using python script (https://github.com/peterthorpe5/intergenic_regions). De novo motif discovery was performed using MEME [80]. The discovered motifs were searched against transcription-factor binding-site databases, such as CollecTF [81], PRODORIC [35], RegTransBase [82], RegPrecise [83], DPINTERACT [84], and Swiss Regulon [85], using Tomtom tool [86]. The listed databases were downloaded from the MEME suite as meme database format, except RegPrecise. For RegPrecise, all transcription factor binding site motifs for Listeriaceae were downloaded and converted to the meme database using sites2meme script from the MEME suite v5.0.5 [80].
Functional enrichment analyses of the DE gene lists GO terms of the genes have been predicted using PANNZER2 [87] with default parameters. The lists of DE genes obtained from the different comparisons and the network gene clusters were tested for the enrichment of GO terms (belonging to the Biological Process ontology) using topGO package [88]. Enrichment for GO terms was tested separately for up-and downregulated genes. The reference set included all the GO term-annotated genes in the genome. The functional categories were considered to be enriched if p-value ≤ 0.05. KAAS (KEGG Automatic Annotation Server) [89] was also used to obtain KO (KEGG Orthology) assignments of genes. ddPCR ddPCR was used to verify RNA-seq results of the barotolerance experiment (n = 18). Three replicate samples of each treatment or strain and their corresponding control samples were always analysed. Expression levels of ten genes: recG, fusA, clpE, hly, agrB, ftsE, mscL, p A, dnaK, and murA, were quanti ed from ScottA samples treated with 200 MPa and recovered for 24 h, or treated with 400 MPa and recovered for 10 min. Expression levels of seven genes: recG, fusA, clpE, hly, agrB, ftsE, and mscL were quanti ed using ScottA samples treated with 400 MPa and recovered for 24 h. Primers (Table S15) were designed using Primer3Plus [90] and manufactured by Integrated DNA Technologies. The protocol included gDNA removal and RT-PCR steps as described previously [91]. To be able to compare expression levels of different samples, expression of the target genes (cDNA copies/µl) was normalized using concentrations of two stably expressed genes: recG and fusA. To help comparison to RNA-seq, the results were expressed as log 2 (gene concentration in treated sample/gene concentration in control sample) values.
Deletion of lmo1013 and lmo2229 genes from L. monocytogenes EGDe genome The bacterial strains and plasmids used in the present study are listed in Table 1. Culture media utilized for the cultivation of E. coli and L. monocytogenes were Luria-Bertani (LB) broth (Sigma Aldrich, St. Louis, Missouri, United States) and brain heart infusion (BHI) broth (Oxoid), respectively. E. coli EC10B chemically competent cells were prepared with the CaCl 2 method [92] and L. monocytogenes EGDe electrocompetent cells were obtained as described earlier [93]. The selective antibiotics and their concentrations were used as follows: kanamycin, 50 µg/mL for E. coli; erythromycin, 250 µg/mL for E. coli and 5 µg/mL for L. monocytogenes; chloramphenicol, 7.5 µg/ml for L. monocytogenes. Two L. monocytogenes EGDe mutants were generated by chromosomal mutagenesis of lmo1013 and lmo2229 genes, based on the system composed of pORI280AD recombinant vector and pVE6007 helper plasmid, following the protocol provided by Monk et al. [93].
The oligonucleotide primers used to amplify the anking regions that shared 20 bp overlapping ends with PstI linearized pORI280 vector (Thermo Scienti c) and 20 bp overlapping ends between them are presented in Table 2. After gel extraction and puri cation (FavorPrep™ GEL/PCR Puri cation Kit, Favorgen), the gene deletion anking fragments were ligated into the cut pORI280 backbone by Gibson assembly (2X Gibson Assembly® Ultra Master Mix, Synthetic Genomics Inc) reaction (pORI280AD). Following chemically competent E. coli EC10B cells' transformation with the recombinant vector by heat shock, colony PCR was employed to screen for transformants. This was done by AD fragment ampli cation from the total DNA released from the heat treated cells (94 °C, 15 min.) with KAPA Taq DNA Polymerase (Sigma Aldrich) and primers pairs lmo1013_AB_fwd / lmo1013_CD_rev and lmo2229_AB_fwd / lmo2229_CD_rev.
Further, the electroporation of L. monocytogenes EGDe and the site-directed mutagenesis of lmo1013 and lmo2229 genes were performed according to the protocol of Monk et al. [93]. Gene deletion screening was performed by Colony PCR using the primers pairs AB_fwd/CD_rev to amplify the genomic region that encompasses both anking fragments and gene of interest.
Gene deletion from L. monocytogenes EGDe genome was further con rmed by DNA sequencing. The same methods were used for genome sequencing and assembly of wild-type EGDe, Δlmo1013, and Δlmo2229 as described previously [12]. In addition, the reads were mapped to reference L. monocytogenes EGDe genome and using Bowtie2 v2.3.4.3 [75] and visualized using IGV v2.4.19 [97] to focus on deletion regions.
Resistance of L. monocytogenes EGDe mutants to HPP The L. monocytogenes EGDe strains were grown in TSB with 0.6% yeast extract (TSB-YE) to early stationary phase, prepared as described previously in High pressure processing for RNA-seq method and then subjected to HPP at 300 and 350 MPa, 8  PA, AIN, DB, TL, CUR, and NB conceived and designed the study. TMR, LG, TL, FIB, LGG, BN, AIN, and DB collected the samples and performed the pressure treatments and viable cell count. FIB, PC, and LGG performed gene deletion experiments. ICD, MA, and PL performed the bioinformatics analyses. AY performed RNA extraction, sequencing preparation, and ddPCR. LP organized NGS assays. PC extracted DNA. ICD, MA, and FIB drafted the manuscript. All authors have read, commented, and approved the nal manuscript.