Strand specific RNA-sequencing and membrane lipid profiling reveals growth phase-dependent cold stress response mechanisms in Listeria monocytogenes

The human pathogen Listeria monocytogenes continues to pose a challenge in the food industry, where it is known to contaminate ready-to-eat foods and grow during refrigerated storage. Increased knowledge of the cold-stress response of this pathogen will enhance the ability to control it in the food-supply-chain. This study utilized strand-specific RNA sequencing and whole cell fatty acid (FA) profiling to characterize the bacterium’s cold stress response. RNA and FAs were extracted from a cold-tolerant strain at five time points between early lag phase and late stationary-phase, both at 4°C and 20°C. Overall, more genes (1.3×) were suppressed than induced at 4°C. Late stationary-phase cells exhibited the greatest number (n = 1,431) and magnitude (>1,000-fold) of differentially expressed genes (>2-fold, p<0.05) in response to cold. A core set of 22 genes was upregulated at all growth phases, including nine genes required for branched-chain fatty acid (BCFA) synthesis, the osmolyte transporter genes opuCBCD, and the internalin A and D genes. Genes suppressed at 4°C were largely associated with cobalamin (B12) biosynthesis or the production/export of cell wall components. Antisense transcription accounted for up to 1.6% of total mapped reads with higher levels (2.5×) observed at 4°C than 20°C. The greatest number of upregulated antisense transcripts at 4°C occurred in early lag phase, however, at both temperatures, antisense expression levels were highest in late stationary-phase cells. Cold-induced FA membrane changes included a 15% increase in the proportion of BCFAs and a 15% transient increase in unsaturated FAs between lag and exponential phase. These increases probably reduced the membrane phase transition temperature until optimal levels of BCFAs could be produced. Collectively, this research provides new information regarding cold-induced membrane composition changes in L. monocytogenes, the growth-phase dependency of its cold-stress regulon, and the active roles of antisense transcripts in regulating its cold stress response.


Introduction
The human pathogen Listeria monocytogenes continues to pose a challenge in the food industry, where it is known to contaminate ready-to-eat foods and to grow during refrigerated storage. Ingestion of L. monocytogenes by susceptible individuals can cause potentially fatal food-borne infections with mortality rates as high as 40% reported [1]. The ubiquitous nature of this pathogen makes it difficult to eliminate from our food systems however, post-processing levels of L. monocytogenes contamination are often low [2][3][4] and unlikely to cause disease [5,6]. Furthermore, the risks associated with L. monocytogenes contamination can be addressed by heating foods to an appropriate temperature [7]. Therefore, the ability of L. monocytogenes to grow to unsafe levels in ready-to-eat foods during the shelf-life of products, pose the greatest concern to consumer health.
When subjected to low temperatures, all living cells will experience similar challenges that stem from a reduction in molecular dynamics, which leads to decreased rates of diffusion and structural changes to molecular structures [8]. Cold-stress adaptation mechanisms among microbes include membrane compositional changes, compatible solute uptake, and the synthesis of nucleic acid stabilization proteins, DNA-unwinding enzymes, and general stress response and cold shock proteins [9,10]. Currently, most of our knowledge regarding bacterial cold stress response (CSR) mechanisms comes from the model organisms Escherichia coli and Bacillus subtilis. However, unlike L. monocytogenes, neither of these bacteria can multiply at temperatures close to 0˚C. In the last decade, L. monocytogenes outbreaks associated with ready-to-eat foods in both North America and Europe have increased interest in the specific mechanisms employed by this psychrotrophic pathogen to adapt and grow at low temperatures.
The CSR of L. monocytogenes has been studied using microarrays [11][12][13], quantitative real-time PCR (qPCR) [14][15][16][17], mutagenesis [18][19][20][21][22][23][24] and other genetic and proteomic techniques [25][26][27]. Collectively, these methods have identified a large pool of genes with putative or known roles in cold tolerance. However, although a small portion of qPCR-based studies have looked at the differential expression (DE) of select genes during the early stages of the L. monocytogenes CSR (i.e., lag phase) [16,17], global transcriptome studies have only been conducted on cold-adapted exponential-and stationary-phase cells to date, leaving much to be discovered regarding the initial CSR of L. monocytogenes. Furthermore, since CSRs become more pronounced in response to more dramatic changes in temperature, many studies have focused on abrupt temperature down shifts from 37˚C to 15˚C or lower [11,14,15,17,28]. Such conditions may not represent the food industry, where more commonly L. monocytogenes is typically transferred from an ambient-temperature environment to a food product or plant environment that is maintained at 4-10˚C. Though exposure to ambient temperature (20-25˚C) can also be considered a low temperature stress [15,29], the maximum growth rate (μ max ) of L. monocytogenes' is approximately 10× faster at 20˚C than at 4˚C whereas there is much less of a difference between 15˚C and 37˚C (~4.5×) [30]. Moreover, at 37˚C, the transcriptional landscape of L. monocytogenes undergoes a drastic change in preparation for intracellular survival [31]. physicochemical cues. These transcripts are often referred to as riboswitches or thermosensors. One of the best-known examples in L. monocytogenes is the riboswitch that blocks the translation of the major virulence regulator prfA at low temperatures (<30˚C) [33]. Trans-encoded sRNAs, on the other hand, are not located adjacent to their target and share only limited complementarity allowing them to regulate multiple mRNAs. Lastly, asRNAs are transcribed from the DNA strand opposite of a gene and thus have perfect complementarity. asRNAs can be short (<100 nt) or long (>1000 nt) and in some cases correspond to the 5'-or 3'-extension of an mRNA transcribed from an adjacent gene [34]. The prevalence of genes with reported asRNA ranges from 20-75% depending on the organism [35,36]. To date, more than 100 asR-NAs have been described in L. monocytogenes [31,[36][37][38]. Compared with the numbers reported for other bacteria this number seems rather low. However, it will presumably rise with further study.
The objective of this study was to gain a comprehensive view of the L. monocytogenes CSR, with an emphasis on elucidating mRNA and asRNA expression patterns across multiple growth phases following cold stress. To gain a better understanding of the mechanisms employed by L. monocytogenes to grow during refrigerated storage, we aimed to characterize the CSR from cells during the initial lag phase following cold stress, throughout exponential growth and long-term stability at low temperatures, using conditions similar to those present in a food-contamination scenario. To accomplish these goals, strand-specific RNA sequencing and membrane lipid profiling were conducted on L. monocytogenes cell cultures at five distinct growth phases, following a downshift in temperature from 20˚C to 4˚C. As the phenotypic and genetic properties of L. monocytogenes isolates differ, we chose to elucidate the CSR mechanisms of a strain with relevance to the food industry and consumer health: a fast cold-growing, serotype 1/2a, food plant isolate containing full length versions of important virulence genes: inlA, inlB, inlC, prfA, plcA, hly, mpl, actA, plcB. Here, we report novel findings regarding coldinduced membrane compositional changes in L. monocytogenes, the growth phase-dependent cold-stress regulon, and the active roles of antisense transcripts in regulating the CSR of this pathogen.

Culture conditions and time point selection
A previously evaluated L. monocytogenes environmental isolate (BioSample SAMN05256775) from a food-processing plant in British Columbia (Lm1, serotype 1/2a, sequence type 7) was selected for use in this study [39]. The strain displayed enhanced cold tolerance (μ max at 4˚C) relative to a large collection of isolates [39]. Bacterial cultures were grown for 24 h in brain heart infusion broth (BHIB; Difco, Fisher Scientific, Canada) at 20˚C, re-suspended in pretempered BHIB at a cell density of 10 7 CFU/ml and incubated at 4 or 20˚C. RNA and lipids were extracted at five time points from cells grown at 4˚C (treatment, T) and 20˚C (control, C). Three biological replicates were conducted for each treatment and control sample. Each time point corresponded to a specific growth phase (G) (Fig 1). The time points were as follows: G1 -early lag phase, G2 -transition to exponential growth phase, G3 -mid-exponential growth phase, G4 -transition to stationary phase, and G5 -late-stationary phase. G1 corresponded to 20% of the lag phase duration at each temperature. The timing was determined by modelling previously obtained growth curve data using the Baranyi and Roberts model [40] in DMfit (v3.5) (http://browser.combase.cc/DMFit.aspx) [30]. G2 was marked by a doubling in cell numbers, confirmed using plate counts. G3 corresponded to mid-exponential growth (10 8 CFU/ml), and G4 corresponded to the transition to stationary phase; both time points were identified spectrophotometrically (A 600nm ) and confirmed with plate counts. Lastly, G5 was selected to correspond to G4 plus 8× the lag phase duration (h) at each temperature.

Fatty acid analysis
At each of the five growth phases for cells grown at both 20˚C and 4˚C, cultures were pelleted (10-45 mg net weight), rinsed twice with one ml phosphate buffered saline (PBS, Fisher Scientific), and stored at -80˚C. One sample per time point was collected. The frozen pellets were later sent to MIDI labs (Microbial ID, Inc., Newark DE, USA) where the cell lipids were extracted, followed by methylation of the fatty acids (FAs), and then loaded onto a gas chromatograph for analysis. Fatty acid methyl ester (FAME) profiles were then generated and analyzed using Sherlock 1 pattern recognition software.

RNA isolation and sequencing
At each sample point, growth in cultures was stopped by adding 10% phenol:chloroform (Fisher Scientific) in ethanol solution pre-chilled to -80˚C in a 1:10 volume to the sample. The tubes were vortexed briefly and then centrifuged immediately for 10 min at 4,696×g and 0˚C. The supernatants were removed and the resulting pellets were stored at -80˚C for less than 2 weeks. Total RNA was isolated and purified using the PowerMicrobiome™ RNA Isolation kit (MO BIO Laboratories, CA, USA) per the manufacturer's protocol. RNA integrity numbers (RINs) were determined using the 2100 Bioanalyzer (Agilent, CA, USA). Samples with a RIN between 9.7 and 10 and an RNA concentration >100 ng/μl were sent to Genome Québec (Montréal, QC, Canada) for rRNA-depleted Illumina TruSeq RNA library prep and TruSeq stranded total RNA 100 bp paired-end sequencing on the Illumina HiSeq 2000 platform.

RNA-seq data analysis
Sequencing quality was assessed using FastQC [41], and Illumina adapter sequences and lowquality base pairs were removed using default parameters in Trimmomatic v 0.36 [42]. After the removal of low-quality reads, 22-39 million reads remained for each dataset. Reads were mapped to the complete sequenced genome of L. monocytogenes EGD (NCBI RefSeq NC_022568.1) using Bowtie 2 v 2.2.6 [43] and allowing zero mismatches. L. monocytogenes EGD was selected as the reference genome as it allowed for the highest percent of successfully mapped reads compared to the more commonly employed reference strains EGDe and 10403S. The mapping efficiency ranged from 97.8-99.6% for individual reads and 92.1-95.6% for successfully paired reads. BAM alignment files were used as input for read counting using featureCounts v 1.5.0-p1 [44]. The default counting mode 'union' was used, and separate count files were generated for sense and antisense transcripts. Differential expression (DE) analyses were performed using DESeq2 [45] in R v 3.1.1 [46], and the DE was reported as log 2 fold changes. p-values were adjusted by the DESeq2 default Benjamini-Hochberg (BH) adjustment method and genes with a >2-fold (>1 log 2 ) change in expression and an adjusted p-value < 0.05 were considered as DE. A principle component analysis (PCA) was conducted using DESeq2 to determine the overall reproducibility of our RNA-seq biological replicates. To evaluate the overall transcription levels for individual genes at both 4˚C and 20˚C, raw read counts were normalized based on library size. No p-values were calculated for genes with low mean normalized counts, zero counts, or extreme outlier counts.

Clustering of gene expression profiles
Clustering of genes with similar expression patterns across the five growth phases, was performed using the Mfuzz package [47] in R. Default filtering and standardization methods (based on standard deviation) were applied to the log 2 data from all five growth phases. Soft clustering was then performed using the fuzzy c-means algorithm and the following parameters: c = 20, m = 1.25. Genes were assigned to clusters of given expression patterns across the growth-phases based on having a membership value >0.5.

Functional categorization of differentially expressed genes
To investigate the roles of DE genes at each growth phase, we determined the overrepresentation of molecular and biological pathways (Fisher exact test, p<0.05) using SmartTables based on the BioCyc database (https://biocyc.org/) [48]. SmartTables were also used to map L. monocytogenes EGD genes to the equivalent genes in L. monocytogenes 10403S for which the BioCyc database contains transcription regulation and gene ontology information. Additional enrichment analyses were then performed (Fisher's exact test, p<0.05) for gene ontology (GO) terms, and genes regulated by certain transcription regulators (i.e., σ B , PrfA, RpoD, VirR, MogR, CodY, HrcA and CtsR).
Quantitative PCR validation of RNA-seq data RNA-seq results were validated using qPCR amplification of three genes: 1) leuA which exhibited >4-fold higher expression at 4˚C at all five growth phases, 2) cspB which exhibited >4-fold lower expression at 4˚C at all growth phases, and 3) recJ which was chosen as a reference gene because it showed very little variation in expression across all growth phases at either temperature (S2 Table). Up to 1 μg of RNA from each of our 30 samples was reverse transcribed using the QuantiTect Reverse Transcription Kit (Qiagen, Valencia, CA) per the manufacturer's protocol. Primers were designed using Primer3 Plus (Table 1) and our draft whole  [31], who reported that Listeria spp. express more than 98% of their ORFs at 30˚C and 37˚C. Three of the 30 sequenced samples were excluded from further analysis as they were deemed to be outliers as visualized in the PCA plot in Fig 2. The excluded samples included one biological replicate each from the T4, C4, and C5 treatments. Two of the three samples were from G4 which represents the transition from exponential growth to stationary phase. Given the short window of time available for extracting RNA from this very specific physiological phase, it is likely that these samples were processed either too early or too late relative to the other two biological replicates. The excluded C5 sample belonged to the same biological replicate as the excluded C4 sample. The PCA plot (Fig 2) also shows that the level of transcriptome variance is much greater (67%) among growth phases than between the two temperatures (9%), highlighting the importance of growth phase selection in determining the outcomes of an experiment.
Overall, 1.3× more genes were suppressed than induced in L. monocytogenes under cold stress (Fig 3). This finding contrasts with the results reported in a microarray study conducted by Chan et al. [11], where 30-40% fewer genes were downregulated than upregulated in L. monocytogenes at 4˚C compared to 37˚C. However, in agreement with their work, we observed that cold-adapted stationary-phase cells exhibited the greatest number of DE genes (G5, Fig 3). From here on, the DE of genes at 4˚C vs. 20˚C will be discussed with respect to the growth phases (G1-G5) and C1-C5 and T1-T5 will be reserved for when referring to overall levels of transcription or changes in the membrane lipid profiles that occur at 20˚C and 4˚C, respectively. Fig 4 showed the pairs of growth phases that shared the highest and lowest numbers of the same up or downregulated genes at 4˚C relative to 20˚C. G4 had the lowest number of upregulated genes in common with the other growth phases, whereas G1 and G2 (G1:G2) and G1:G5 shared many of the same upregulated genes. With respect to downregulated genes, G4 and G5 had the smallest in common with G2 and G3 whereas G1:G2 and G1:G5 shared a high number of downregulated genes.
To confirm our DE findings, we used RT-qPCR to quantify the expression levels of a gene that was upregulated for >4 fold in response to cold at all growth phases (cspB), and a gene that was downregulated for >4 fold in response to cold at all growth phases (leuA), according to the RNA-seq data. We observed a strong positive correlation (R 2 = 0.96, y = 1.20x-0.41) between the mRNA levels detected by the two methods. Nevertheless, the fold changes calculated using RNA-seq data were consistently larger than those obtained using RT-qPCR (S1 Fig).

Core set of genes induced by cold stress
The DE analyses revealed a core set of 22 genes induced at 4˚C relative to 20˚C at all five growth phases (Table 2). These genes fall into the category of cold acclimation proteins (CAPs), which are defined as proteins encoded by genes that show continuous increased expression throughout prolonged cold-stress exposure [8,50]. Among these CAPs were nine genes involved in the production of branched-chain amino acids (BCAAs) (ilvDBHC-leuABCD-ilvA). These BCAA synthesis genes were statistically overrepresented among the upregulated genes at all growth phases ( Table 3). The BCAAs leucine, isoleucine, and valine, are critical for BCFA production, with isoleucine being the precursor for anteiso BCFAs. Listeria and similar bacteria, such as Bacillus species, alter their membrane lipid composition in response to temperature downshifts, to include higher percentages of branched-chain fatty acids (BCFAs) with anteiso configurations [24,[51][52][53]. Anteiso BCFAs have lower melting points than their iso-branched and straight-chain counterparts. They are thus more effective at increasing membrane fluidity, which is needed to maintain membrane-transport functions at low temperatures.
Other genes with increased expression in response to cold stress at all growth phases, included those encoding the L-carnitine transporter OpuC (opuCB-opuCC-opuCD), internalins A and D, a ribosomal protein (RpmF), a histidine biosynthesis protein (HisE), a GNAT family acetyltransferase (LMON_0625), and Veg protein (LMON_0187). The roles of carnitine  and ribosomal proteins in bacterial adaptation to certain stresses have been previously reported [18,[54][55][56]. In B. subtilis Veg protein has been shown to activate extracellular matrix genes and contribute to biofilm formation [57]. In L. monocytogenes, increased expression of Veg has been observed following exposure to acid and cold stress [11,58], but its function remains unknown. In addition to the upregulation of inlA and inlD at all growth phases at 4˚C, inlH and two genes encoding internalin-like proteins (LMON_2407, LMON_0611) were upregulated at four growth phases ( Table 2). The most recognized roles of internalins A and B are in L. monocytogenes' virulence. Nevertheless, the Listeria spp. family of internalin proteins all share a leucinerich-repeat domain that allows them to bind structurally unrelated ligands, thereby implicating them in a wide range of functions [59]. Recent studies have shown that isolates containing a full-length version of inlA exhibit enhanced cold tolerance relative to those with a truncated version [39,60]. Further research is needed to confirm the potential role(s) of internalins in the L. monocytogenes CSR.
Many more genes (n = 104) were upregulated at four growth phases at 4˚C relative to 20˚C, with the majority being induced at G1, G2, G3, and G5 (G1-G3:G5, n = 70). These genes shared similar DE patterns across all growth phases (Fig 5. Cluster 5 and 12), where levels of DE dramatically decrease at G4. Present in this group were genes involved in arginine biosynthesis (dapE, argGHJ), general stress response (LMON_0515, LMON_2234, LMON_2696, LMON_2770), FA synthesis (LMON_0351), and phosphotransferase system (PTS) uptake of cellobiose (LMON_2800-2803) and mannose (LMON_0785-0788) ( Table 2). Additionally, the genes encoding succinate-semialdehyde dehydrogenase (LMON_0920), 11 genes from the pentose phosphate pathway (LMON_0350-0360), and two sets of genes (LMON_0354-0356, LMON_2718-2720) encoding dihydroxyacetone kinase (DhaKLM) were upregulated in these growth phases (Table 2). An additional 100 genes displayed induced transcription at three growth phases, with the largest number co-upregulated at G1-G2:G5 (n = 29), followed by G1-G3 (n = 24). The similarities between these sets of growth phases are visible in Fig 4. At G1-G2:G5, L. monocytogenes growth is inhibited, and such conditions would be expected to induce genes encode proteins associated with 1) secondary membrane-transport systems that are less energy demanding, 2) the utilization of alternative energy sources, and 3) cell detoxification. It appears that the demand for these proteins is further increased in cold stressed cells as we observed the upregulation of genes involved in utilizing ethanolamine as an alternative energy source (LMON_1167-1180), oligopeptides transport (oppC), transcription regulation (5 genes), general stress response (clpC), and ribosomal subunit assembly (LMON_2523). At G1-G3, several of the shared upregulated genes encoded reductases (LMON_1496, LMON_2403, LMON_0798, LMON_0643, LMON_2843) with known or putative roles in maintaining redox balance within the cell and managing oxidative stress.
Genes involved in the biosynthesis and export of peptidoglycan, teichoic acids, cell wall proteins, and membrane lipids were also suppressed at 4˚C (Table 4). Both peptidoglycan and teichoic acid biosynthesis stem from the precursor molecule UDP-N-acetyl-α-D-glucosamine (UDP-GlcNAc) and require the attachment of D-alanine residues, and for peptidoglycan also glutamate. Correspondingly, genes involved in UDP-GlcNAc biosynthesis and transport, alanine attachment (dltABCD), and glutamate transport (LMON_0849-0850) and attachment (racE) were downregulated. Other suppressed genes included several fatty-acid-coA ligases, lipid export proteins, and enoyl-acyl-carrier proteins and reductases, which collectively assist in the production and export of phospholipids, cell wall-associated lipoproteins, and lipoteichoic acids. These findings suggest a reduced rate of peptidoglycan/cell envelope turnover at 4˚C relative to 20˚C. This reduced turnover rate may reflect the reduced cellular growth rate at this temperature, or conservation of carbohydrates for alternative uses.
Other genes suppressed at 4˚C had roles in localization and pathogenicity, heme biosynthesis, tRNA processing and charging, nucleotide degradation and salvage, and carbohydrate transport and catabolism (Table 4). Not surprising, many of these genes also have alternative roles in cell wall biogenesis. Notably, over 50 genes from PTS operons were downregulated at one or more growth phases, far exceeding the number of PTS genes upregulated in response to cold stress. PTSs utilize phosphate to facilitate the uptake of simple sugars and thus consume more energy than other membrane kinases with the same sugar specificity [69]. This suggests that L. monocytogenes may benefit from employing alternative sugar uptake systems when exposed to cold stress to conserve energy for more critical CSR mechanisms.

L. monocytogenes cold-stress response at individual growth phases
Early lag phase (G1). Previous studies have described the lag phase following cold stress as the acclimatization phase. In this phase, bacteria suppress bulk protein synthesis, and dramatically increase production of transient cold induced stress response proteins (CIPs) [8]. Once the cells are cold adapted, CIP synthesis ends and bulk protein synthesis and cell growth resume. Depending on the degree of cold stress, the types and expression levels of CIPs will differ. Some of the functions associated with CIPs include DNA unwinding, nucleic acid stabilization, and compatible solute uptake [70,71]. A total of 96 genes were uniquely induced at G1, accounting for 25% of the genes upregulated at this growth phase. Most of these genes had DE patterns that belonged to clusters 14, 16, or 20 (Fig 5), and a large number encoded transcription regulators (n = 19), including HrcA, CtsR, DegU, and YycF. Both HrcA and CtsR are negative regulators of several heatshock and general stress-response genes, a number of which were also upregulated at G1 (grpE, dnaK, clpB, gmpA, mscAB, clpC, clpE, and clpP). Many of these upregulated genes, including groEL and groES, encode chaperone proteases capable of degrading the misfolded or damaged proteins that would probably occur following a rapid temperature downshift. By contrast, the transcription regulator DegU contributes to motility, growth at high temperatures, and efficient biofilm formation [72]. Co-transcribed with degU is yviA, which encodes a DegV family protein with putative functions in FA transport or metabolism [72,73]. As we did not observe upregulation of genes known to be activated by DegU, the role of this operon in early cold adaptation may be primarily associated with the functions of YviA. In support of this hypothesis, YycF, a transcription regulator of FA biosynthesis and cell wall metabolism genes [74], was also upregulated at G1.
Twenty-nine genes were solely induced at G1 and G2, highlighting their importance in the earlier CSR stages of L. monocytogenes. Included among these genes were those encoding the low temperature requirement protein LtrC, cold-shock protein CspL, regulatory protein RecX, and several proteins associated with phosphate transport (PstCABB-PhoU-PstS).
Forty-seven genes were upregulated at both G1 and G5 (late-stationary phase), many of which followed the DE pattern shown in cluster 16 ( Fig 5). As discussed, these growth phases both represent stages in cold stress survival, where the population size is static and their induced genes are probably specific to surviving adverse conditions. Among these genes were those encoding heat-shock proteins (clpP, clpE, groEL, LMON_2710-2711), DNA gyrase (gyrA), an SOS-response protein (lexA), an RNA helicase (secA), and seven transcription regulators. Interestingly, although these genes were upregulated in response to cold stress in this study, several of them were also reported by van der Veen et al. [75] as induced in L. monocytogenes following 3 min of heat stress (48˚C). As both stresses are caused by a rapid change temperature which causes molecular structures to take on different conformations; it makes sense they would share similarly upregulated genes with functions in DNA and RNA repair (gyrA, lexA, secA), and protein degradation (clpP, clpE) or stabilization (groEL).
Transition to exponential growth phase (G2). At T2, cells have had an adequate amount of time to repair cold-induced damage to molecular structures, and to reprogram their cellular machinery to support growth at low temperatures. As T2 shares characteristics with both the lag and exponential phases, only 22 genes were uniquely upregulated at G2, and only 21 genes were uniquely downregulated. The exclusively induced genes encoded the iron binding protein, Fri; a cold-shock protein, CspD; and a DNA repair and SOS response protein, RecA. Most of these genes exhibited slightly less than a 2-fold increase in expression at G1, maximum increased expression at G2, and no change in expression or decreased expression from G3-G5, as is characteristic of cluster 10 ( Fig 5). The upregulation of fri, cspD, and recA have been previously reported in L. monocytogenes in response to cold stress [16,26,27,50]. Among the genes upregulated at G2, there was a significant (p<0.05) overrepresentation of those involved in glycerol metabolism and the synthesis of arginine, histidine, and BCAAs (Table 3). Among the most induced genes (>36 fold) were those encoding the OpuC carnitine uptake system (opuCABCD), hypothetical proteins, and succinate-semialdehyde dehydrogenase (LMON_0920) which was also strongly induced at G1 (Table 5).
Mid-exponential growth phase (G3). A total of 47 genes were uniquely induced during exponential growth at 4˚C, with most genes exhibiting the DE pattern shown in cluster 13 ( Fig 5). Aside from a few genes with roles in iron transport (fhuC, fhuB, LMON_2440-2441) and oxidoreductase activity (LMON_0560, LMON_2031, pdhD), most are not known to share any common functions. Among G3-upregulated genes, there was an overrepresentation (p<0.05) of genes involved in glycerol metabolism, aerobic respiration, iron transport, and the biosynthesis of arginine, histidine, and BCAAs (Table 3). Other upregulated genes included those involved in the pentose phosphate pathway (LMON_0349-0352, LMON_0647), osmolyte uptake (opuCABCD), and 12 genes fell within PTS operons that facilitate the uptake of mannose, cellobiose, and beta-glucoside. Genes highly induced at G3 (>16-fold) were associated with iron (efeUOB) and heme transport (isdCAEF), as well as protein export (tatAC) ( Table 5). As expected, many of the proteins previously induced during cold acclimatization (G1-G2) were downregulated during exponential growth.
Transition to stationary phase (G4). Upon entry into stationary phase, cold-adapted cells downregulated almost 3× more genes than they upregulated (Fig 3). Unique to G4 was the upregulation of 53 genes with roles in glycine betaine transport (gbuAB), tRNA processing * indicates <0.50 cluster membership; G1-G5 refer to differential expression in L. monocytogenes cells grown at 4˚C relative to 20˚C and across five specific growth phases (see Fig 1); Blue shading indicates genes with significantly increased gene expression (> 1 log 2 , p<0.05) and yellow shading indicates genes with significantly decreased expression (< -1 log 2 , p<0.05) at 4˚C relative to 20˚C; Bolded values indicate significant (p<0.05) differential expression changes. Genes shaded in yellow are highly expressed at more than one growth phase. https://doi.org/10.1371/journal.pone.0180123.t005 (LMON_T39, T60, T66), nucleotide biosynthesis (pyrP, pyrB, pyrAa, pyrR, LMON_1838), and maintaining a rod cell shape (mreBD), among others. Most of these genes displayed DE patterns that belonged to clusters 8 or 19 (Fig 5), in which DE peaks at G4 and falls dramatically at G5. The five most strongly induced genes (>16-fold) at this growth phase encoded three isoleucine biosynthesis proteins (ilvDBH) and two cell wall-associated proteins (LMON_2534, LMON_1312) ( Table 5). G4 and G5 also shared a set of exclusively upregulated genes (n = 23) with functions in nucleotide biosynthesis (purAHNFE, pyrDII, carB, pyrC), tRNA processing (LMON_T51), and cell-wall recycling (LMON_0029, LMON_2776), among others. All genes exhibited patterns of DE that belonged to cluster 1 or 15 (Fig 5). Nucleotide biosynthesis genes are likely activated at these growth phases to accommodate the increased demand for ribonucleic acids imposed by the large abundance of strongly upregulated genes at G5.
Late-stationary phase (G5). Cold-adapted stationary-phase cells exhibited both the largest number of DE expressed genes and the largest magnitude of expression changes (Fig 3). Similar results were noted by Chan et al. [11] in their L. monocytogenes microarray study, which compared the cold-stress regulons of exponential and stationary-phase cells. More than 50% of G5-upregulated genes exhibited >4-fold increased expression and 17 genes, predominantly associated with PTS mediated sugar uptake systems, were upregulated between 100and >1,000-fold. Many of the genes strongly induced by cold stress at this growth phase, also exhibited some of the highest expression levels in stationary-phase cells at 20˚C, highlighting their overall importance during this growth phase and their enhanced importance in coldadapted stationary-phase cells. Additionally, 50-60% of G5 genes were exclusively DE at this growth phase. The exclusively induced genes had roles in tryptophan (trpBCDFG, aroE) and ATP synthesis (atpACDFGH, adk, LMON_0091), creatine and rhamnose degradation, and ribosome and phage-related processes. Other genes encoded motility-associated proteins (motB, cheV, flaA), an osmosensitive K+ channel histidine kinase (kdpCDE), an Na(+) H(+) antiporter (LMON_2393-2396), the transcription repressor CodY, and low temperature requirement protein B (ltrB). Many genes associated with transcription and translation were also strongly induced at G5 (polC, dnaC, dnaI, topB, rbfA, mfd), including 28 transcription regulators, and 14 genes involved in tRNA processing (S1 Table).

L. monocytogenes cold stress regulon: A comparison with previous studies
Numerous previous studies have elucidated a large pool of genes with putative or known roles in the L. monocytogenes CSR [11-13, 16, 20, 27]. In this section, we will determine whether these genes are expressed in a cold-tolerant strain, whether a significant induction occurs when the downshift is from 20˚C to 4˚C (rather than from 37˚C, as in most prior studies), and whether any growth-phase dependencies exist.
Osmolyte and oligopeptide uptake. In both L. monocytogenes and B. subtilis, low temperature and high osmolarity stress induce the intracellular accumulation of solutes and short peptides. These solutes and peptides function as osmoprotectants that facilitate cell growth under stressful conditions [18,26,54,[76][77][78][79]. In L. monocytogenes, carnitine and glycine betaine are the predominant solutes that accumulate [78][79][80], and their transport systems are encoded by the opuCABCD and gbuABC operons, respectively. Chan et al. [11] reported induction of opuCABCD and gbuC in cold-adapted L. monocytogenes exponential-but not stationary-phase cells, while Durack et al. [12] did not find that either system was induced in late exponential early stationary-phase cells. In the present study, opuCBCD was upregulated in all growth phases, while opuCA was upregulated at G1-G4 (Table 6). Regardless, all genes in the OpuC operon belonged to the DE pattern shown in cluster 10 (Fig 5), exhibiting the highest level of induction at G2 (up to 48-fold) and the lowest level at G5, which agrees with previous work [11,12]. Similar to the findings of Chan et al. [11] gbuA and gbuB were only upregulated at G4 (Table 6).
Oligopeptide uptake in L. monocytogenes is mediated by the membrane permease OppA (oppABCDF) which transports peptides containing up to eight residues [23]. This transporter appears to contribute to L. monocytogenes cold tolerance as an oppA null mutant displayed reduced growth at 5˚C in BHI medium [23]. Previously, Durack et al. [12] observed upregulation of oppBCF in late exponential early stationary-phase cells at 4˚C while Chan et al. [11] found that only oppA levels were elevated in stationary-phase cells. Contrary to these reports, we did not observe any notable DE of the opp operon, except for a 2.6-fold induction of oppF at G5 (Table 6). While the operon may be necessary for low temperature growth, it does not appear to have increased transcription at 4˚C relative to 20˚C in our cold tolerant strain.
RNA and DNA repair. At low temperatures, both DNA replication and transcription are hindered by cold-induced changes in nucleic acid structures. L. monocytogenes appears to respond to these challenges by upregulating genes encoding topoisomerases and DNA gyrases, which help maintain the superhelical tension of DNA; RNA helicases, which unwind secondary RNA structures; and exo-and endonucleases, which function in nucleic acid repair [11,12,15]. Markkula et al. [15] reported the upregulation of four DEAD-box RNA helicase genes (lmo0866, lmo1246, lmo1450, lmo1722) in L. monocytogenes for up to 7 h following a downshift from 37˚C to 5˚C. They also found that, compared with the wildtype strain, null mutants of lmo0866, lmo1450, and lmo1722 displayed restricted growth at 3˚C. Similarly, we observed upregulation of lmo1450 (LMON_1513) and lmo1722 (LMON_1787) at G1-G4, and lmo0866 (LMON_0869) at G1:G4 (Table 6). Also in agreement with Markkula et al. [15], lmo1246 (LMON_1305) appeared to be the least important of the four genes, with <2-fold induced expression at G1-G3.
Several DNA gyrase and topoisomerase genes were also upregulated in our study. Notably, topA was upregulated at G4, and topB and parC were upregulated at G5 (Table 6). In L. monocytogenes, the DNA gyrase genes gyrB and gyrA are in a four-gene operon, along with a gene encoding a hypothetical protein, and recF, a DNA recombination and repair protein. All four genes were upregulated at G1, and gyrA was also induced at G5. Upregulation of another DNA recombination and repair protein, recA, was also observed at G2.
Among the nucleic acid repair-proteins induced in this study were ribonuclease J1 (LMON_1037) at G1, an endoribonuclease (LMON_0846) at G2-G3, two nuclease subunits (LMON_2342-2343) at G5, and an excinuclease (LMON_0848) at G4-G5 (Table 6). Additionally, an X family DNA polymerase (LMON_1225) was induced at G1:G5. These polymerases synthesize unusual DNA structures that may serve as indicators of the induction of DNA repair [81][82][83]. A DNA-binding protein (LMON_2812) with a putative role in oxidativedamage protection [84] was also upregulated at G2 and G5. Lastly, a hypothetical protein (LMON_1317) containing a nudix hydrolase domain was upregulated at G1-G3:G5 (Table 6). Some members of this protein family are known to degrade oxidatively damaged nucleoside di-and triphosphates, while other members control the levels of metabolic intermediates and signaling compounds [85].
Overall, genes associated with RNA and DNA repair were predominantly cold-induced in lag-and stationary-phase cells (G1 and G5). This is expected given that cold-induced damage to nucleic acids is highest directly following cold stress, and nutrient-depleted stationary-phase cultures experience higher mutation rates and increased levels of oxidative damage [86][87][88][89]. Regulatory elements. Many mutant characterization and transcriptome studies have evaluated the roles of alternative sigma factors (σ B , σ C , σ H , σ L ), two-component regulatory systems (TCRSs), and negative regulators in the L. monocytogenes CSR. Of the four alternative sigma factors, σ B has been the most extensively studied. It positively regulates over 100 genes when the organism enters stationary phase or is subjected to environmental stresses including low pH, high salt, or carbon starvation [90][91][92][93][94]. Although induced expression of sigB has been observed in L. monocytogenes up to 12 h following cold stress [11,17,19], the σ B regulon is predominantly downregulated during cold stress [11,12], demonstrating that increased sigB expression does not necessarily correlate with increased σ B activity. Furthermore, compared with wildtype strains, sigB null mutants show little to no difference in cold tolerance [16,19,20,95], indicating that sigB is not critical for cold-stress survival. In the present study, the DE pattern of sigB followed that of cluster 10 (Fig 5), with a maximum 1.5-fold increase in expression at G2 (Table 6). Similar expression patterns were observed for three additional sigB operon genes (rsbWVX) while the first four genes (rsbRSTU) of the eight-gene operon were not induced under cold stress. Despite low levels of sigB induction, we found that σ B -dependent genes were overrepresented (p<0.05) among the cold-induced genes at all growth phases ( Table 7). Examples of such genes include fri, bsh, dapE, uspA, gabD, opuCABCD, mpoBACD, phoU, csbD, ltrC, and hrcA. Although some of these genes are solely activated by σ B , others are co-regulated by other unknown or known transcription factors such as σ L , σ H , and PrfA [16,20,[96][97][98]. As previously mentioned, σ B is also activated upon entry into stationary phase [99][100][101]. Consistent with this fact, we observed that sigB transcript levels during growth at 20˚C increased from 4.5k mapped reads at C1 to 15k at C5 (S2 Table).
The alternative sigma factor σ L has been shown to contribute to the ability of L. monocytogenes to tolerate cold [20,102], osmotic, and acid stress [102,103]. This sigma factor positively regulates >400 genes, including those involved in cell envelope synthesis, motility, and PTS sugar uptake and catabolism [28,104]. Although sigL induction has been reported in exponential and late exponential early stationary-phase cells of cold-adapted L. monocytogenes [12,27,102], deleting sigL does not impact the ability of L. monocytogenes ability to grow at cold temperatures [20,28]. In the present study, the DE pattern of sigL followed that of cluster 5 (Fig 5), with baseline expression at G1-G4 and <2-fold increased expression at G5 (Table 6). Sixteen genes previously identified as being positively regulated by σ L in L. monocytogenes at 3˚C [28] also exhibited similar DE patterns. The remaining alternative sigma factors, σ H and σ C , contribute to the survival of L. monocytogenes under alkaline and heat stress, respectively [105][106][107]. However, these sigma factors have limited roles in cold tolerance [20]. Consistent with previous findings, we did not observe the upregulation of sigH or sigC in our study.
RpoD is the principal RNA polymerase sigma factor in L. monocytogenes and regulates housekeeping genes associated with ribosome structure, protein synthesis, and rRNA and tRNA [108]. Like the DE pattern for sigL, the DE pattern of rpoD followed that of cluster 5 (Fig 5), with a maximum 1.5-fold increase in expression at G5 (Table 6). Accordingly, genes positively regulated by RpoD were overrepresented (p<0.05) among the downregulated genes at all growth phases except G4. Many genes in the RpoD regulon are co-activated by PrfA, and correspondingly, PrfA-regulated genes were also overrepresented among the downregulated genes at G1-G3:G5 (Table 7). Genes co-regulated by RpoD and PrfA include prfA, plcAB, hly, mpl, inlC, and actA, all of which are associated with virulence [109]. Other studies have also reported lower levels of PrfA-regulated virulence genes in cold-adapted exponential and stationary-phase cells [11,12,110]. In the present study, the DE pattern of prfA fit that of cluster 17 (Fig 5), with 2-fold and 40-fold increased expression at G3 and G5, respectively (Table 6). Interestingly, although the transcription of prfA increased dramatically at G5, the expression levels of many PrfA-dependent virulence genes remained strongly downregulated. Researchers have proposed that PrfA activity, but not prfA transcription, is inhibited by unphosphorylated forms of PTS permeases that occur during active sugar transport; these PTS permeases may bind directly to PrfA [111][112][113]. Our results support this hypothesis, as many PTS operons and other carbohydrate uptake and catabolism genes were upregulated at G5 and belonged to the same cluster profile as prfA (Fig 5, cluster 17). Whether PrfA is subsequently involved in regulating these genes is still unclear.
CodY is a pleiotropic transcriptional regulator that actively represses the transcription of genes involved in amino acid metabolism, nitrogen assimilation, mobility and chemotaxis, and sugar uptake, among others [114]. In agreement with previous work [11], we observed the highest induction (2.4-fold) of codY in at G5 following the DE pattern shown in cluster 3 ( Fig 5).
DegU is another well-known response regulator that regulates the expression of motility-, virulence-, and biofilm-related genes in L. monocytogenes [29,72,115]. Our results showed that degU induction was greatest (3-fold) at G1 (Table 6), with a DE pattern consistent with that of cluster 16 (Fig 5). This probably explains why increased expression of degU was not observed in previous cold stress transcriptome studies of L. monocytogenes, which focused on exponential-and stationary-phase cells [11,12]. Despite the induction of degU at G1, the flagella operon (LMON_0680-0694) that the DegU protein regulates remained suppressed, likely due to co-induction of the gene encoding the transcriptional repressor MogR, which shared the same DE pattern as degU. In B. subtilis, DegU is part of the two-component system DegS/U, which regulates the expression of genes encoding various extracellular enzymes [116]. This suggests that DegU may also contribute to the regulation of extracellular enzymes in L. monocytogenes that facilitate cold growth.
The negative regulators HrcA and CtsR repress expression of several genes for heat-shock proteins and cellular protein quality control (e.g. dnaK, grpE, groES, groEL, clpB, clpC, clpE, clpP) in L. monocytogenes and other bacteria [117][118][119]. Increased expression of many of these genes has also been reported in response to salt, cold, and ethanol stress [27,120,121]. In the present study, we found that hrcA and ctsR, like degU and mogR, were both maximally upregulated at G1, as were many of the genes they regulate. This was somewhat surprising given the inverse relationship expected between the DE patterns of transcription repressors and their gene targets. However, CtsR and HrcA regulons are known to have a considerable amount of overlap with PrfA and σ B regulons, demonstrating the complexity and fine-tuning abilities of bacterial regulatory networks, which allow bacteria to survive a wide-range of conditions [98,122,123].
Two-component signaling systems (TCSs) are also important regulators of bacterial stress responses and typically consist of a transmembrane sensor histidine kinase (HK), and a cognate cytoplasmic response regulator [124][125][126][127]. The sequenced genome of L. monocytogenes EGD-e contains 16 known TCSs [128]. Using mutant characterization, Pöntinen et al. [129] showed that only the HKs LisK and YycG of the LisKR and YycGF TCSs, respectively, are important for L. monocytogenes growth under cold stress. However, increased expression of lisK was not observed in either the Pöntinen et al. [129] study or the present study, highlighting that the importance of a gene does not necessarily correlate with its level of induction. HKs with increased expression in our study included cesK and LMON_1166 at G1, and kdpD, LMON_1166, and virS at G5 (Table 6). Of these, kdpD, which encodes the HK for an osmosensitive K+ channel, exhibited the largest increase in expression (18.6-fold). KdpD, together with its response regulator KdpE, controls the expression of a high-affinity K+ translocating ATPase [130,131]. Among the many K+ transport systems in bacteria, Kdp has the highest affinity for K+ and it is only expressed when other systems are unable to meet the cell's K+ needs. Our results therefore suggest that activation of Kdp is an important mechanism used by L. monocytogenes to survive long-term exposure to cold stress.
Ribosome functions. When bacteria are subjected to a temperature downshift, their ribosome structures become compromised, resulting in translation inhibition and a prolonged lag phase. In fact, inhibition of ribosomal functions induces CSR proteins [132]. In E. coli, three ribosome-associated proteins (IF2, CsdA, RbfA) are required for protein synthesis and subsequent cell growth at low temperatures [133][134][135]. In contrast, 22 other ribosomal proteins are not believed to be essential for the growth of E. coli or B. subtilis at low temperatures [136,137]. In L. monocytogenes, the role of ribosomes and their associated proteins in cold stress remains largely unknown. In a study by Durack et al. [12], ribosome protein genes were strongly activated in osmo-and cold-adapted L. monocytogenes cells. In our study, of the 58 ribosomal proteins identified in L. monocytogenes EGD, 22 were upregulated exclusively at G5 in the cold tolerant-strain we studied. Among these genes were those encoding L. monocytogenes homologs of IF2 (infB) and RbfA, previously mentioned for their roles in stabilizing cold-sensitive ribosomes in E. coli. An additional two ribosomal protein genes, ctc and rpmF, were significantly upregulated at all growth phases, and a ribosomal subunit interface protein (LMON_2523) was induced at G1-G2:G5 (Table 6). In B. subtilis, ctc is induced in response to osmotic, heat, and oxidative stress [138,139], and ithas similarly been linked to osmo-and cold-tolerance in L. monocytogenes [12,56].
In bacteria, rRNA and ribosomal protein synthesis are tightly controlled, to meet the translational needs of the cell. The increased transcription of ribosome proteins in cold-adapted stationary-phase cells may reflect an increased demand for protein synthesis, as suggested by the large number of strongly upregulated genes at G5. Alternatively, ribosomal proteins have been shown to participate in extra-ribosomal functions, as independent polypeptides with roles in transcription and DNA repair [140][141][142]. Thus, they might contribute to the L. monocytogenes CSR in yet undetermined ways.
Cold-stress proteins. CSPs are a conserved family of small (~70 aa) proteins containing a nucleic acid-binding domain. Found in many prokaryotic and eukaryotic organisms, CSPs can act as transcriptional activators, antiterminators, or as RNA chaperones that enhance translation at low temperatures by blocking the development of secondary mRNA structures [10,70,71]. Three CSPs have been identified in L. monocytogenes and are listed here in the order of functional importance: CspL>CspD>CspB [143]. Furthermore, L. monocytogenes CSPs appear to be only induced during the early stages following cold stress [11,12,143]. In the present study, cspB was downregulated at all growth phases, reaching a maximum 228-fold decrease at G5, while cspD and cspL were upregulated at G1-G2 and then exhibited no change or decreased expression at the later growth phases (Table 6).
Additional proteins with putative roles in the L. monocytogenes cold-stress response. A few other genes have also been associated with the L. monocytogenes CSR. Zheng and Kathariou [144,145] identified three low temperature requirement proteins (LtrA, LtrB, and LtrC) necessary for growth at cold temperatures. However, varying results have been reported regarding the expression of these genes at low temperatures. Pieta et al. [14] observed higher levels of ltrC transcripts at 7˚C than at 37˚C in exponential-phase cells, whereas other researchers saw either no difference in expression or decreased expression in exponential-and stationary-phase cells of cold-adapted L. monocytogenes [11,12,27]. Increased expression of ltrC has also been reported in L. monocytogenes EGD-e, directly following an upshift in temperature from 37˚C to 48˚C [146]. Importantly, we observed increased expression of ltrC only in the early stages following cold stress (G1 and G2), which may partly explain why no differences were seen in studies that analyzed exponential-and stationary-phase cells. As for other low temperature requirement protein genes, ltrB was upregulated at G5 while no notable changes were seen for ltrA.
The fri (flp) gene, which encodes ferritin, is also commonly discussed in cold stress studies of L. monocytogenes and is hypothesized to play a role in iron storage [96,146,147]. Previous studies have reported that fri is induced upon entry into stationary phase in L. monocytogenes subjected to low temperatures [27,50,96]. Additionally, fri null mutants display reduced growth at 4˚C in BHIB and increased sensitivity to oxidative stress [146,148]. In our study, fri was exclusively induced at G2 (Table 6). However, in agreement with the findings of Polidoro et al. [96], the number of fri transcripts at 20˚C was nine times higher in stationary-phase cells than in lag-phase cells.
Liu et al. [22] identified a membrane-associated phosphohydrolase (pgpH) as the interrupted gene in an L. monocytogenes cold-sensitive transposon mutant. Compared with the parent strain, this mutant also showed increased intracellular levels of the phosphorylated guanosine nucleotide (p)ppGpp, suggesting that PgpH may be critical for (p)ppGpp degradation at low temperatures. Cellular (p)ppGpp inhibits RNA synthesis when a shortage of amino acids is present. Such conditions can occur during cold stress, because of structural damage to membrane transporters and intracellular enzymes. Arguedas-Villa et al. [17] reported increased pgpH expression in cold-tolerant but not cold-sensitive strains of L. monocytogenes at 4˚C compared to 37˚C. Though we also evaluated the gene expression of a cold-tolerant strain, we observed no significant induction of pgpH (LMON_1529) at 4˚C relative to 20˚C.
The induction of flagella biosynthesis and motility-related genes has frequently been reported in L. monocytogenes following a temperature downshift [11,14,27]. However, the reliability of these observations is debatable as 37˚C was uniformly used as the control temperature and L. monocytogenes is not typically motile above 30˚C [149,150]. Recently, Cordero et al. [13] reported that fast cold-growing L. monocytogenes strains are less motile than slow cold-growing strains. They hypothesized that low motility might allow cold-tolerant strains to proliferate more rapidly at low temperature. In the present study, our cold-tolerant strain exhibited decreased expression of most flagella operons at all growth phases. However, five motility-specific genes (LMON_0691-0695) were upregulated at G5, including motB, gmaR, cheV, and flaA ( Table 6). The role of motility genes in prolonged cold-stress survival remain speculative, but one hypothesis is that it is beneficial for cells to be able to move to environments more suitable in terms of nutrition or other requirements.

Cold-induced membrane lipid composition changes
To gain a better understanding of both the timing associated with cold-induced membrane FA changes in L. monocytogenes and the types of changes that occur, we analyzed FAs extracted from cells at the same time points used in our RNA-seq experiment.
To date, the membrane FA profile of L. monocytogenes under cold stresses has primarily been investigated using late exponential to early stationary-phase cells [24,51,151,152]. In most cases, researchers have observed a~20% increase in the relative proportion of a-C15:0 among cold-grown cells (5-10˚C), with maximum levels ranging from 66-80%. Correspondingly, a-C17:0 levels decrease by 20-25% with minimum levels ranging from 3-14%.
In the present study, a-C15:0 levels increased from 46% at T1 to a maximum of 70% at T4 (Fig 6A), whereas levels remained constant around 50% in 20˚C-grown cells. Thus, in this study, the first one to look at cold-induced membrane FA changes at multiple growth phases, we show that L. monocytogenes continues to make alterations until it transitions into stationary phase (T4), at which point minimal further adjustments are made. Correspondingly, a-C17:0 levels decreased from 12% at T1 to 4% at T4, whereas at 20˚C they ranged from 13-18% (Fig 6C). Levels of i-C16:0 and iC15:0 also decreased by 2.3% and 9%, respectively at 4˚C (Fig 6B). Additional changes in the BCFA composition of cold-stressed cells included a 2.5-2.9% increase in the proportion of i-C14:0 and the appearance of i-C13:0 and a-C13:0.
As previously discussed, anteiso BCFAs are synthesized from the BCAAs isoleucine, and iso BCFAs are synthesized from leucine and valine [52]. These BCAAs are produced by proteins encoded by a nine-gene operon (LMON_2054-2062). At 4˚C, we found that all genes in this operon were induced >4-fold at all growth phases. Interestingly, genes from the BKD operon (LMON_1432-1437) which encodes enzymes that convert BCAAs into BCFAs, showed either no change in expression or decreased expression at 4˚C. Other researchers have also noted a lack of bkd induction in L. monocytogenes under cold-stress conditions [24,51]. A BKD-independent pathway may exist, in which BCFA synthesis uses branched-chain α-keto acids instead of branched-chain acyl-CoA esters as primers [52]. Alternatively, Nickel and colleagues [153] showed that in B. subtilis the promoter regions of five operons, including BKD, are not cold-inducible but exhibited increased mRNA stability at low temperatures.
Shortening of fatty acid chain lengths. Following a temperature downshift, L. monocytogenes incorporates FAs with shorter chain lengths to further decrease the phase-transition temperature of its membrane [156]. Overall, we observed an 8% increase in C14 FAs at 4˚C relative to 20˚C, which became noticeable at T3 (Fig 7A). Moreover, by T5, 92% of Lm1's membrane FAs contained fewer than 15 carbons; in 20˚C grown cells, this percentage was 75%. Our results suggest that L. monocytogenes begins degrading unfavorable membrane phospholipids during lag phase but cannot synthesize its preferred FAs until growth resumes. This poses the question of how the cell retains adequate membrane fluidity during the early stages of cold stress. The answer to this question appears to involve the desaturation of existing FAs.
Increase in unsaturated fatty acids. Several bacterial species including E. coli, some species belonging to the phylum cyanobacteria, and some species belonging to the genus Bacillus, increase their proportion of membrane unsaturated FAs (UFAs) directly following a temperature downshift [157][158][159]. This process occurs rapidly due to the presence of membraneassociated desaturases that introduce double bonds into preexisting saturated FAs (SFAs). For example, in E. coli, cis-vaccenate is produced 30 seconds following a shift from 42˚C to 24˚C [160]. To date, the roles of UFAs in the L. monocytogenes CSR have been largely overshadowed by interests in BCFA modulation. Although some papers have reported the absence of UFAs in the L. monocytogenes membrane [24,[161][162][163], others have detected their presence in small amounts under various conditions [51,151,152,162,[164][165][166][167][168]. In a study of L. innocua, levels of C18:1 FAs increased to 6.2% following a downshift from 35˚C to 10˚C [169]. Similarly, another study found that the Bacillus megaterium membrane contained a maximum of 22% C16:1 cis-Δ 9 (palmitoleic acid) at 10˚C, whereas the maximum was 2% at 35˚C [53]. In the present study, the proportions of C16:1 cis-Δ 9 and C18:1 cis-Δ 9 (oleic acid) peaked at T3 (7%) and T2 (20%), respectively (Fig 7B). By contrast, no difference was detected between the proportions of oleic and palmitoleic acid in T4 and C4 cells, or T5 and C5 cells. Very small percentages (<1.2%) of other UFAs were also detected in our samples, but no notable differences were observed between the temperatures treatments (S3 Table). The fact that membrane UFA proportions only increased in cold-stressed lag-and exponential-phase cells likely explains why studies that have investigated late exponential-and stationary-phase cells have failed to detect any changes.
Unlike the rapid conversion of SFAs to UFAs, switching from iso to anteiso BCFAs and from longer to shorter acyl chain FAs requires de novo synthesis of whole lipid molecules by cytoplasmic enzymes that are usually linked to growth [170]. Sato and Murata [171] showed that following a 10-15˚C downshift in temperature, cyanobacterial cells could only resume growth and FA biosynthesis once a certain degree of membrane unsaturation was reached. Based on our findings, we suggest that in response to cold stress, L. monocytogenes converts C18:0 and C16:0 to oleic and palmitoleic acid, respectively, to rapidly lower its membrane phase-transition temperature. Fig 8 shows how UFA levels appear to compensate for a-C15:0 until its optimal levels are reached at T4. The centre placement of the cis-double bond in oleic and palmitoleic acid makes these FAs highly efficient at increasing membrane fluidity. For example, the average melting points of oleic and palmitoleic acid are 13˚C and 1.22˚C, respectively, compared to 24˚C and 52˚C for a-C15:0 and i-C15:0, respectively [172]. While oleic and palmitoleic acids are effective at rapidly increasing the membrane fluidity of L. monocytogenes, this membrane profile is probably less stable than one that includes a high proportion of a-C15:0. In support of this idea, Juneja and Davidson [173] showed that when provided with an external source of oleic acid, L. monocytogenes was able to increase the relative proportion of this acid in its membrane from 0.72 to 28%. However, these cells were highly susceptible to salt and several antimicrobials.
In addition to increasing in response to cold, UFA levels also increased up to 11% in lag and exponential-phase cells grown at 20˚C (Fig 7B), suggesting a role for UFAs in regular cell growth. In Saccharomyces cerevisiae, oleic acid and ergosterol supplementation mitigate oxidative stress [174]. Similarly, bacteria transferred from a stationary-phase environment into a fresh oxygenated medium also experience oxidative stress [175] and may benefit from increased levels of membrane UFAs.
In Bacillus, phospholipid desaturases, encoded by des, are induced upon cold stress and controlled by a two-component regulator that senses changes in membrane fluidity [176,177]. These desaturases introduce double bonds at the Δ 5 or Δ 10 positions of acyl chains attached to existing phospholipids. Pseudomonas has a similar desaturase system that can introduce double bonds at the Δ 9 position [178]. Currently, no lipid desaturases have been identified in Listeria. However, Vadyvaloo et al. [166] showed that when L. monocytogenes is treated with a desaturase inhibitor, it exhibits decreased levels of membrane UFAs, demonstrating that such a system does exist. The synthesis of oleic (18:1 cis-Δ 9 ) and palmitoletic (16:1 cis-Δ 9 ) acids by L. monocytogenes suggests that it contains a desaturate system like that present in Pseudomonas. In the present study, we demonstrate the presence of additional UFAs, which contain double bonds in several other locations (S3 Table) highlighting the possibility that more than one desaturate system exists in the bacterium.

Comparison of antisense transcription at 20˚C and 4˚C
A total of 92-261k paired-end (PE) reads were successfully mapped to the antisense strand of L. monocytogenes EGD ORFs, accounting for 0.61-1.63% of the total number of mapped reads. Overall, 70% of genes had >10 mapped PE asRNA reads, however, only 56 genes had >1,000 asRNA mapped PE reads with a maximum of 40,000 PE reads observed for a LMON_R5 (vs. a maximum of 591,000 observed for a single mRNA transcript). These results are in line with previous findings from other bacteria, from which asRNAs were detected for 20-75% of genes [35,[179][180][181].
Overall, more asRNA transcripts were upregulated than downregulated under cold stress (Fig 3). The largest number of upregulated asRNA transcripts was observed at G1 (Fig 3), which makes sense given that asRNAs act as gene regulators. At T1, cells are actively restructuring their transcriptional activity and cellular processes to prepare for growth in their new environment [175]. Around 10% of the asRNAs currently described for L. monocytogenes are thought to be involved in regulating transcription regulators [34]. This idea is supported by the overall low levels of antisense expression observed in this study.
While many asRNA transcripts were upregulated at G1, total antisense transcription was generally highest in late stationary-phase cells at both temperatures (Table 8), highlighting the importance of antisense regulation at this growth phase. Compared with regulatory proteins, asRNA transcripts offer the advantage of being able to provide a rapid connection between quorum sensing and direct destabilization of target mRNAs; moreover, as they can act post- transcriptionally, they can more tightly control the repression of proteins under environmental stress conditions [34,182,183]. Finally, unnecessary regulatory RNAs can be cleared quickly and via a process that requires less energy consumption than the removal of regulatory proteins. This may in part explain the abundance of asRNAs during late stationary-phase, in which energy sources are limited. Most genes with no or <10 asRNA transcripts encoded tRNAs (S2 Table). By contrast, rRNA had the highest levels of antisense transcription (Table 8). Similarly, Wehner et al. [38] have reported high levels of antisense transcription for rRNA in L. monocytogenes under intracellular growth conditions. Given the critical importance of ribosomes in cell functioning, it makes sense that bacteria would take advantage of the tight means of regulation offered by antisense transcripts. Other highly expressed antisense transcripts targeted genes encoding the cell-shape determining protein, MreBH; flagellar biosynthesis protein, FliP; magnesium and cobalt transport protein, CorA; and two transcription regulators (LMON_0384, LMON_0738), among others (Table 8). Again, most of these genes exhibited maximum levels of asRNA transcription in stationary-phase cells with higher levels evident at 4˚C than 20˚C.
Some L. monocytogenes genes exhibited high levels of antisense transcription with no or very low levels of mRNA transcription. Not surprisingly, this occurred when an antisense transcript was an extension of a 3' or 5' untranslated region (UTR) of an adjacent gene (Fig 9A, 9B and 9F), suggesting that the asRNA probably prevented RNA polymerase from binding to the promoter region of the gene on the opposite strand. On the contrary, antisense transcripts that appeared to be individually transcribed, frequently had high levels of mRNA expression ( Fig  9E, 9G and 9M). These forms of asRNA have been shown to alter transcript stability by forming a sense/antisense RNA duplex leading to RNase-mediated degradation [183][184][185][186] or by stabilizing mRNA transcripts by inducing cleavage of unstable polycistronic transcripts [187][188][189]. Future research aimed at validating the mechanisms and functions of specific L. monocytogenes antisense transcripts will further enhance our understanding of gene regulation in this pathogen and possibly lead to the development of novel intervention strategies that can be used in the food industry.

Conclusions
Here we present results from the first time-course study to investigate the transcriptional response and associated cold-induced membrane FA changes of L. monocytogenes, from early growth phases through late stationary-phase. To the best of our knowledge this is also the first study to look at asRNA expression in L. monocytogenes during cold stress and across multiple growth phases. Our results revealed that the L. monocytogenes transcriptomic response to cold stress is most active during late stationary-phase survival. Similarly, we show that the L. monocytogenes cold-stress regulon differs greatly across growth phases, highlighting the importance of carefully selecting appropriate time points when designing and conducting transcriptome studies.
Overall, more genes were suppressed than induced in L. monocytogenes under cold stress conditions. A core set of 22 genes was upregulated at all growth phases, including nine genes required for BCFA synthesis, the osmolyte transporter genes opuCBCD, and genes encoding the internalins A and D. This reflects the cell's need to synthesize BCFAs to maintain membrane fluidity at low temperatures, and to support proper protein-folding through the uptake of structural-stabilizing osmolytes. Genes suppressed at 4˚C were largely associated with cobalamin (B12) biosynthesis or the production/export of cell wall components. While is unclear how L. monocytogenes may benefit from downregulating cobalamin biosynthesis genes during cold stress, a reduced rate of peptidoglycan/cell envelope turnover at 4˚C relative to 20˚C may reflect the reduced cellular growth rate at this temperature, or allow carbohydrates to be conserved for alternative uses.
Notable cold-induced membrane FA changes included a 15% increase in the proportion of BCFAs and a 15% transient increase in UFAs between the lag and exponential phases. Such information may be useful for improving intervention strategies that target the food industry, as increased membrane UFA levels are known to increase the susceptibility of L. monocytogenes to salt and several antimicrobials.
Overall, around 70% of L. monocytogenes genes exhibited antisense transcription but only a small portion of antisense transcripts exhibited expression levels comparable to that of mRNA.  Table 6). The coverage maps shown are from stationary phase cultures grown at 4˚C. Panels A-F shows individual paired-end, reads while panel G shows overall coverage trends for rRNA (5, 16, and 23S). Reads with the same start and end alignment positions are overlaid and depicted in green. Areas shaded in blue and orange denote sense and antisense transcripts of target genes, respectively. https://doi.org/10.1371/journal.pone.0180123.g009 On average, antisense transcription was higher at 4˚C than at 20˚C, highlighting the importance of antisense regulation in the L. monocytogenes CSR. The largest number of upregulated antisense transcripts was observed during early lag phase; however, at both temperatures antisense transcription was generally highest in late stationary-phase cells. Stationary-phase cells likely benefit from the tight control of protein expression that is offered by antisense RNA, thereby reducing the overall energy requirement of the cell while faced with nutrient-limiting conditions.
Collectively, our results reveal novel gene expression patterns and membrane FA alterations that occur in L. monocytogenes following cold stress, and highlight the abundance of antisense transcription in this microorganism at both 20˚C and 4˚C. We believe that this research will serve as a platform for future cold-stress studies by facilitating the selection of candidate genes and antisense transcripts for further functional validation.
Supporting information S1 Table. Differential expression (log 2 ) and associated cluster memberships of L. monocytogenes EGD ORFs and antisense RNA at 4˚C vs 20˚C. Genes belonging to a single cluster exhibited !0.50 membership. For ORFs with 0.2-0.5 membership to multiple clusters, the clusters are listed in decreasing order of membership. A Ã denotes ORFs with <0.5 membership to one cluster while a "0" denotes ORFs with <0.2 membership to one or more clusters. NC = average normalized count for control treatment. (XLSX) S2 Table. Normalized paired-end read counts for all L. monocytogenes EGD ORFs and antisense RNAs at 4 and 20˚C. Paired-end read counts were averaged across biological replicates. "as" denotes antisense RNA transcripts. (XLSX) S3 Table. Fatty acid membrane profiles of L. monocytogenes at 4˚C and 20˚C. Values are presented as percentages relative to the total amount of extracted fatty acids. (XLSX) S1 Fig. Correlation between differential expression (log2) levels obtained using RNAsequencing and quantitative real-time PCR (qPCR). Points represent the average differential expression values for the genes cspB (▲) and leuA (•) at all five growth phases evaluated in this study. The y-axis represents the differential expression levels obtained using RNA sequencing while the x-axis represents the levels obtained using qPCR.