Transcriptomic time-series analysis of cold- and heat-shock response in psychrotrophic lactic acid bacteria

doi:10.21203/rs.3.rs-34513/v1

Download PDF

Research article

Transcriptomic time-series analysis of cold- and heat-shock response in psychrotrophic lactic acid bacteria

https://doi.org/10.21203/rs.3.rs-34513/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 07 Jan, 2021

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background:

Psychrotrophic lactic acid bacteria (LAB) species are the dominant species in microbiota of cold-stored modified-atmosphere-packaged food products and they are the main cause of food spoilage. But still, the cold- and heat-shock response of the spoilage-related psychrotrophic lactic acid bacteria has not been studied. Here, to study cold- and heat-shock response of spoilage lactic acid bacteria, we performed time-series RNA-seq for Le. gelidum, Lc. piscium and P. oligofermentans using temperatures of 0 °C, 4 °C, 14 °C, 25 °C and 28 °C.

Results:

We showed that the cold-shock protein A (cspA) gene was the main cold-shock protein gene among cold-shock protein genes in all three species. Our results indicated DEAD-box RNA helicase genes (cshA, cshB) play a critical role in cold-shock response in psychrotrophic LAB. In addition, several RNase genes were also involved in cold-shock response in Lc. piscium and P. oligofermentans. Moreover, gene network inference analysis provided candidate genes involved in cold-shock response. Ribosomal proteins, tRNA modification, rRNA modification, and ABC and efflux MFS transporter genes clustered with cold-shock response genes in all three species, which was a strong indication that these genes would be part of cold-shock response machinery. Heat-shock treatment caused upregulation of Clp protease and chaperone genes in all three species and we were able to identify transcription binding site motifs for heat-shock response genes in Le. gelidum and Lc. piscium. Finally, we showed that food spoilage-related genes were upregulated at cold temperatures.

Conclusions:

The results of this study provide new insights into a better understanding of the cold- and heat-shock response in psychrotrophic LAB. In addition, candidate genes involved in cold- and heat-shock response predicted using gene network inference analysis could be used as a target for future studies.

Epigenetics & Genomics

Food spoilage

RNA-seq

gene network inference

time-series

differential gene expression

gene interaction

stress

cold tolerant

Lactic acid bacteria (LAB) are a group of gram-positive bacteria with a wide range of phenotypic and genomic features [1]. On one hand, LAB communities play an important role in fermented foods during the production stage and they can be also used as food preservatives [2]. On the other hand, psychrotrophic LAB cause food spoilage in cold-stored modified-atmosphere-packaged (MAP) food products, since they are able to prevail in the MAP food environment [3]. LAB species composition and their relative abundance depend on the nature of food products and preservation technology [4, 5]. However, two LAB species, Leuconostoc gelidum and Lactococcus piscium, have frequently been predominating at the end of shelf life in a variety of packaged and refrigerated foods of animal and plant origin [6, 7, 8, 9]. Spoilage communities also contain less abundant and slower growing species, such as Paucilactobacillus oligofermentans (former Lactobacillus oligofermentans), the role of which in food spoilage is unclear [10, 11]. We have been working on these three different LAB for several years; we have sequenced the genomes of these species [12, 13, 14, 15] and also analysed their gene expression patterns in growth experiments [14, 15, 16]. Since with these species reverse genetics methods are very slow (for Le. gelidum) or not working at all so far (for Lc. piscium and P. oligofermentans), detailed omics analysis is the only efficient way to study them. Understanding gene expression mechanisms of these spoilage LAB is critical, since MAP technology with combined cold storage has increased its popularity for preservation of minimally processed fresh foods, and better understanding of LAB’s genomics and especially mechanisms of cold-shock and stress adaptation is a key to discover new methods to control spoilage activities.

There are three main categories of bacteria based on their ability of growth at different temperatures; thermophiles, mesophiles, and psychrophiles, they are able to grow at high, intermediate, and low temperatures, respectively [17, 18]. Psychrophiles divide into two: psychrophiles sensu stricto, which optimally grow at 15 °C and psychrotrophic (psychrotolerant), which optimally grow at 20–25 °C [19, 20, 21]. Based on previous published studies, cold-shock protein (CSP), DEAD-box RNA helicase, and ribonuclease (RNase) are commonly known cold-shock response gene families in all three types of bacteria [22, 23, 24, 25]. Similarly, chaperone and Clp family gene families are the common heat-shock response genes in bacteria [18, 26]. To our knowledge, cold- and heat-shock response has been investigated before in mesophilic LAB [27, 28, 29], however, psychrotrophic LAB have not been studied. Here, to investigate both cold- and heat-shock response in spoilage psychrotrophic LAB for the first time, we performed RNA-seq using five temperatures (0 °C, 4 °C, 14 °C, 25 °C and 28 °C) and three time points (5, 35, 185 min) for each temperature. The time points were selected aiming at catching early and also later effects of the shocks, while keeping the sample number still reasonable.

Bacterial growth

Temperature points were selected based on the literature analysis of the biology of psychrotrophic bacteria [19, 20, 21]. The two lowest temperatures used (0 °C and 4 °C) caused a cold-shock effect. To have an additional temperature point between cold-shock and optimum temperature (25 °C), 14 °C was selected. Finally, 28 °C was selected to be heat-shock temperature.

All the bacteria were first grown at 25 °C, and then aliquoted to 5 different temperatures for the certain time (see materials and methods; Figure S1). Le. gelidum and Lc. piscium grew significantly slower at cold-shock temperatures (0 °C and 4 °C) compared to growth in control temperature (25 °C). At 14 °C, notably slower growth was observed only for Le. gelidum indicating Le. gelidum was more sensitive to mild cold-shock temperatures compared to two other species. P. oligofermentans grew slightly slower at cold-shock temperatures (0 °C, 4 °C and 14 °C) compared to growth in control temperature (25 °C), but the difference was not statistically significant. Also, none of the species showed significant growth change at 28 °C compared to control temperature 25 °C (Figure S2).

Leuconostoc gelidum subsp. gasicomitatum LMG 18811

Differentially expressed genes at different temperatures in Le. gelidum

Differential gene expression analysis showed that the number of differentially expressed genes increased by time at different temperatures, except at the mild cold-shock temperature (14 °C), at which the number of differentially expressed genes decreased at the last time point (Figure S3). The general transcriptome view provided by log₂ fold change heatmap of all differentially expressed genes (Figure S4) revealed that upregulation of cold-induced genes increased continuously during time at 0 °C and 4 °C, while continuous increase of expression was not observed at 14 °C. As expected, cold-induced genes were mostly downregulated at 28 °C and on the other hand heat-induced genes were mostly downregulated at 0 °C, 4 °C and 14 °C (Figure S4).

To classify the differentially expressed genes (Table S1), gene ontology (GO) enrichment analysis was performed. The results showed that RNA processing and acetyl-CoA metabolic processes GO terms were enriched for upregulated genes at all cold temperatures (Fig. 1). In addition, enrichment of RNA methylation, and ribosome biogenesis GO terms at cold temperatures suggests that processing of RNA and ribosomal activities take part during cold-shock response. Enriched GO terms for upregulated genes at 4 °C and 14 °C at the 5 minute time points such as peptidoglycan biosynthetic process, cell wall organization, signal transduction by protein phosphorylation, and peptidyl-histidine phosphorylation, indicated that Le. gelidum sensed cold using signal transduction, and cell wall related genes were first overexpressed at cold-shock. At 28 °C, GO terms related to metabolism and protein folding were enriched among upregulated genes (Fig. 1). For downregulated genes, enrichment of ATP synthesis related GO terms was detected at cold temperatures (Figure S5).

Cold-shock, heat-shock and stress related genes in Le. gelidum

We focused on known cold-shock response genes, such as cold-shock proteins, DEAD-box RNA helicases and RNases. Within these genes, cold-shock protein gene cspA (CBL92332.1), DEAD-box ATP-dependent RNA helicase gene cshA (CBL92087.1), RNase R gene rnr (CBL92046.1), and RNase J gene rnj (CBL91141.1) were defined as cold-induced genes in Le. gelidum, since these genes were significantly upregulated at cold temperatures (Fig. 2a,2b,2c). Cold-induced nusA-IF2 operon in E. coli [30] was seen in Le. gelidum, and it (rimP, nusA, ylxR, ribosomal protein L7AE gene, IF-2) was upregulated at cold temperatures. In addition to nusA-IF2 operon, upregulation of translation initiation factor IF-3 (CBL91889.2) was seen at cold temperatures (Fig. 2d). Interestingly, at 0 °C, only cshA was significantly upregulated and the rest of the cold-related genes were not differentially expressed. None of the known cold-shock response genes were upregulated at 14 °C 185 minute time point, although significant upregulation was seen at 35 minute time point (Fig. 2).

Heat-inducible transcription repressor hrcA, chaperone genes (groS, groL, dnaK, and dnaJ), Clp protease genes (clpP, clpE), and chaperone-binding gene grpE were significantly upregulated at the heat-shock temperature (28 °C). The upregulation of most heat-shock genes was not detected at 185 minute time point due to high expression level of heat-shock genes also at 25 °C. The heat-shock genes were significantly downregulated at cold temperatures during the last time point (Fig. 2e).

Depending on the temperature, both upregulated and downregulated stress- related genes were observed for Le. gelidum. At the cold temperatures, GlsB-family stress protein gene (CBL91583.1) and PspC-domain containing stress-responsive transcriptional regulator (CBL90655.1) were significantly upregulated. Interestingly, another GlsB-family stress protein gene (CBL91598.1) was upregulated at 28 °C. Significant downregulation of stress protein genes, such as CsbD-family stress protein genes (CBL91764.1 and CBL91769.1) and universal stress protein gene UspA (CBL90994.1) was detected at cold temperatures (Fig. 2f).

Gene network inference in Le. gelidum

We aimed to identify gene interactions and detect novel cold- and heat-shock response genes. To achieve this, we used a simple guilt-by-association approach by performing gene network inference analysis and gene interaction network based clustering for all differentially expressed genes. According to gene interaction network analysis, 80 clusters that include more than two genes were identified (Table S2). Cold-shock response genes (cspA, cshA, rnr) were seen in three different clusters (Cluster names will be referred as “Cl” from this point): Cl5, Cl22 and Cl34 (Fig. 3, Table S2), from which Cl22 and Cl34 were linked to each other. The genes within Cl5 were significantly enriched for the pseudouridine synthesis GO term. In addition to pseudouridine synthesis genes (CBL92346.1 and CBL91776.1), methyltransferase genes were observed in clusters including cold-related genes indicating that several methyltransferase genes (CBL92090.1, CBL91526.1, CBL90969.1, CBL91168.1, CBL92048.1) were also interacted with cold-shock response genes. The strong interaction of pseudouridine synthesis genes, methyltransferase genes and cold-shock response genes implied that pseudouridine synthesis genes and methyltransferase genes could play a role in cold adaptation. We also observed that two-component system regulatory protein gene yycH (CBL90958.1) was clustered with cshA within Cl5. In addition, two-component system regulatory protein gene yycG (CBL90957.1) and two-component sensor histidine kinase gene hpk4 (CBL92274.1) were located in Cl6, which was linked to Cl5. This indicates these sensor genes might play a role in cold sensing in Le. gelidum. Several transcription factor genes, such as MarR family transcriptional regulator (CBL90793.1), bacteriophage transcriptional regulator (CBL92264.1), HTH-type transcriptional repressor CzrA (CBL92283.1), HTH-type transcriptional regulator GmuR (CBL90950.1), HTH-type transcriptional repressor RghR (CBL91060.1), transcriptional regulator (CBL90737.1), and SsrA-binding protein SmpB (CBL92045.1) were also seen in cold-shock response genes including clusters (Cl5, Cl22 and Cl34) and clusters linked to them.

All heat-shock related genes were clustered in Cl3, and the Cl3 was linked to Cl13 based on gene interaction network (Fig. 3). As expected, the genes within the Cl3 were significantly enriched for protein folding GO term, as most of the heat-shock genes were chaperones. In addition to heat-shock genes, transporter protein genes, such as divalent metal cation transporter MntH (CBL92285.1), folate ECF transporter (CBL92501.1), and chloride transporter (CBL91698.1), were also located in Cl3, which indicates that transportation function was part of heat adaptation. Similarly, Cl13 was significantly enriched for amino acid transmembrane transport GO term. Several transcription regulators were observed within Cl13, specifically carbohydrate diacid transcriptional activator CdaR (CBL91787.1), transcriptional regulator TetR/AcrR family (CBL92525.1), and predicted transcriptional regulator (CBL91595.1), indicating possible role of these factors in the regulation of heat-shock genes.

de novo transcription binding site prediction in Le. gelidum

We wanted to understand whether the genes clustered together by the expression patterns would also be regulated with similar transcription factors. We first checked, whether any known transcription binding site motifs were enriched within all possible upstream regions in the Le. gelidum genome. The result showed that CcpA binding site was the most enriched motif within the whole genome. In addition, we predicted that MalT, galR, GalS, MtrB and rpoD transcription factor binding sites significantly occur in Le. gelidum (Table S3). Since cold-shock protein gene cspA can act as a transcription enhancer by binding to the 5'-ATTGG-3' in the promoter regions of genes [31], we specifically searched for it, and detected 284 upstream regions with the 5'-ATTGG-3' motif (Table S4). Within these 284 regions, there were upstream regions of cold-induced genes such as nusA, beta-lactamase and RNA methyltransferase genes. However, the motif was also present in upstream regions of heat-induced genes such as groL and groS chaperone genes and stress genes (Table S4).

To predict de novo transcription binding sites, motif discovery analysis was done for upstream regions of the upregulated genes for all conditions. Motifs similar to a CcpA binding site were discovered in the upstream regions of upregulated genes at 0 °C, 4 °C and 28 °C (Table S5). At 14 °C, only 35 minute time point upstream regions of upregulated genes showed a statistically significant motif (Table S5), which was not matched with other motifs in transcription factor binding site (TFBS) databases. However, the motif was similar to ribosomal binding site Shine-Dalgarno sequence. No motif was discovered in the upregulated genes at 14 °C 185 minute time point. This was expected since most of the cold-shock response genes were not upregulated at 185 minute time point. In addition, other motifs without database match were discovered in the upstream regions of upregulated genes (Table S5).

To look more closely to co-expressed genes, clusters that were created using gene inference analysis were analysed for the de novo motif discovery. Cold-shock response genes including clusters (Cl5, Cl22, Cl34) were not shown to have any significant motif. Interestingly, a motif with statistically significant p-value was seen in Cl3, which includes heat-shock related genes. Upstream regions of four heat-shock related genes (clpE, groS, hrcA, and clpP) and one hypothetical protein gene were contributing to the construction of the motif (Table S6). Tomtom database match analysis showed that this motif significantly matched with the HrcA motif in RegPrecise database [32]. Upstream regions of Cl4 and Cl15 genes showed similar motifs, which matched with the GalR binding site motif (Table S6). CcpA motif was also significantly matched with the motif, though GalR was the most significant match. This was expected, since several genes in Cl4 and Cl15 were annotated as carbohydrate metabolic process genes (Table S6). In addition, the motifs with significant E-values were discovered in Cl23, Cl35 and Cl42, but no motif database match was seen for discovered motifs (Table S6).

Pathway enrichment and changes of metabolism at different temperatures in Le. gelidum

To observe metabolism changes we performed KEGG pathway enrichment analysis for upregulated and downregulated genes. The results for upregulated genes showed that the two-component system was the only KEGG term enriched at all cold temperatures. Fatty acid biosynthesis and fatty acid metabolism terms were seen up both at 4 and 14 °C. Upregulated genes at 28 °C were mainly enriched for central metabolism KEGG terms, such as glycolysis, starch and sucrose metabolism, and galactose metabolism, which were not observed for the three lower temperatures (Fig. 4). Downregulated genes at 0 and 4 °C were enriched mainly for central metabolism, indicating metabolism was significantly slower at 0 and 4 °C (Figure S6).

Lactococcus piscium MKFS47

Differentially expressed genes at different temperatures in Lc. piscium

Only a couple genes were differentially expressed at cold temperatures at 5 min time point, indicating that gene expression adaptation to cold temperatures happens later in Lc. piscium. However, we observed more than 100 differentially expressed genes at 28 °C already at 5 min time point, which suggests that heat raises a much faster and robust change in gene expression comparing cold-shock treatment. The number of differentially expressed genes increased with time at all temperatures and the number of differentially expressed genes at 0 °C and 4 °C was highest at 185 min, when about half of the genes were differentially expressed (Figure S3, Table S7). Among the species studied, Lc. piscium had the highest number of differentially expressed genes in these conditions. Log₂ fold change heatmap of all differentially expressed genes (Figure S7) showed clear separation between cold- and heat-induced genes. GO enrichment analysis of upregulated genes showed that RNA catabolic process, RNA hydrolysis, nitrogen compound transport, nucleotide salvage, and purine nucleobase metabolic process GO terms were enriched at all cold temperatures, which suggests involvement of RNA catabolism and purine metabolism in cold adaptation. Ribosome and methylation related GO terms were also enriched for upregulated genes in cold temperatures (Fig. 5). Protein folding GO term was the most significantly enriched term for upregulated genes at 28 °C. Glycerol metabolism and stress response GO terms were enriched at 28 °C for upregulated genes (Fig. 5), while these GO terms were enriched for downregulated genes at cold temperatures (Figure S8).

Cold-shock, heat-shock and stress related genes in Lc. piscium

Despite the fact that Lc. piscium is a psychrotrophic bacterium, it harbours only one cold-shock protein gene; cspA (CEN27394.1), which was upregulated at cold temperatures (Fig. 6b). Similarly upregulation of DEAD-box RNA helicase genes (cshb (CEN27700.1), cshA(CEN27648.1)), RNase genes (Ribonuclease HIII, M5, J2, J1, R, Z, Y, P; CEN27237.1, CEN29282.1, CEN28217.1, CEN28898.1, CEN28776.1, CEN28949.1, CEN29011.1, CEN27509.1) and translation initiator factor genes (IF-1 (CEN27470.1), IF-3 (CEN28088.1)) was detected at cold temperatures (Fig. 6a,6c,6d). Cold-induced nusA-IF-2 operon was found in Lc. piscium, but the operon was not upregulated at cold temperatures.

Heat-shock related genes, such as chaperones and heat-inducible transcription repressor hrcA showed significant upregulation at most time points at 28 °C (Fig. 6e), with the simultaneous downregulation of these genes at cold temperatures.

Upregulation of stress genes at cold temperatures was not detected. Moreover, several stress genes, such as universal stress protein gene UspA1 (CEN28155.1), UspA2 (CEN29325.1) and GlsB family stress response protein gene (CEN29488.1), were downregulated at cold temperatures. However, Gls24 family (CEN29584.1) and CsbD family (CEN27327.1) stress-response protein genes were upregulated at 28 °C (Fig. 6f), indicating stress genes were part of the heat-shock reaction.

Gene network inference in Lc. piscium

Based on gene network inference and gene clustering, 87 clusters with more than two genes were predicted (Table S8). Most of the cold-shock response genes (cshA, cshB, Ribonuclease Z, Ribonuclease J) were clustered together in Cl2 (Fig. 7). Ribosomal proteins (50S ribosomal protein L33 and ribosomal protein S1 RpsA) and ribosome biogenesis genes were seen within the Cl2, which indicates that ribosome related genes interacted with cold-shock response genes and might play role in cold adaptation. Interestingly, DNA repair genes, such as recA, recF, and recJ, were clustered together with cold-shock response genes. Similarly to Le. gelidum, methyltransferase and tRNA modification genes (LSU methyltransferase RlmI (CEN28296.1), tRNA uridine 5-carboxymethylaminomethyl modification enzyme MnmG (CEN28086.1), SAM-dependent methyltransferase (CEN28161.1), and tRNA (guanine-N(1)-methyltransferase (CEN28283.1)) interacted with cold-shock response genes. Sensor histidine kinase (CEN29277.1) found in Cl13, which is linked to Cl2, was the only sensor-related gene linked to cold-shock response genes.

All heat-shock genes were clustered together in Cl3 (Fig. 7) and it was not linked with other clusters. Protein folding GO term was the most enriched GO term for the Cl3. TetR family transcriptional regulator gene (CEN27501.1) was located in Cl3 suggesting that it interacts with heat-shock genes.

de novo transcription binding site prediction in Lc. piscium

Within all possible upstream regions in the Lc. piscium genome, rpoD17 motif was the most enriched motif. In addition, CcpA, MtrB, Crp, and CpxR motifs were significantly enriched within the upstream regions of the whole genome (Table S9). We detected 341 upstream regions with the 5'-ATTGG-3' CspA binding site. Upstream regions of both cold-induced genes, such as cshA, cshB, recA, and ribonuclease genes, and heat-induced genes, such as dnaJ and groS, were among those 341 upstream regions (Table S10).

Since Lc. piscium had a large proportion of differentially expressed genes, the motif discovery for the upstream regions of upregulated genes resulted in several significantly enriched motifs. As expected, motif discovery analysis of a large number of upstream regions caused enrichment of the Pribnow box motif (TATAAT). Only two of the discovered motifs were matched with a motif from TFBS database. At 14 °C 35 minute time point, discovered motif matched with PhoP motif from PRODORIC database [33] and at 28 °C 185 minute time point, discovered motif matched with rpoD17 motif from DPInteract database [34] (Table S11). A CtsR-binding site like motif was discovered for upstream regions of upregulated genes at 28 °C 5 minute time point, eventhough de novo motif finding E-value score was not significant. Eleven upstream regions, including heat-shock response genes and transcriptional regulator CtsR, contributed to the construction of the motif (Table S11).

The motif discovery analysis for clusters, which were created using gene inference analysis, showed two significant motifs for Cl2 and Cl3 (cold-shock cluster and heat-shock cluster, respectively) (Table S12). Neither of the discovered motifs were matched with any known transcriptome binding motif. Similarly to upregulated genes at 28 °C at 5 minute time point, a motif matching with a CtsR-binding site like motif, but without significant E-value, was found for Cl3. Thereafter, we merged linked clusters (Cl2, Cl13, Cl34, Cl41, Cl65) for cold-shock response genes and discovered one more motif. All the 53 upstream regions used contributed to the construction of the motif. The discovered motif, however, was not matched with the TFBS database (Table S12).

Pathway enrichment and changes of metabolism at different temperatures in Lc. piscium

KEGG pathway enrichment analysis results did not show any common enriched pathway for upregulated genes at cold temperatures (Fig. 8). Only ribosome KEGG term was enriched at 0 °C and 4 °C, and RNA degradation term at 4 °C and 14 °C. In addition, pathway enrichment analysis using Pathway Tools showed that both cell wall biosynthesis and teichoic acids biosynthesis pathways were enriched for upregulated genes at 0 °C and 4 °C. Glycerololipid metabolism was the only KEGG pathway term enriched for upregulated genes at 28 °C. We also observed central metabolism and amino acid metabolism terms to be enriched for downregulated genes at cold temperatures (Figure S9).

Paucilactobacillus oligofermentans DSM 15707

Differentially expressed genes at different temperatures in P. oligofermentans

P. oligofermentans had a few differentially expressed genes (Table S13) at the first time point at cold-shock temperatures. However, differential expression of 42 genes was detected at the first time point at heat-shock temperature 28 °C, which suggested that gene expression reaction for this species was faster at heat-shock conditions compared to cold temperatures. The number of differentially expressed genes was increased during time at 0 °C and 4 °C, while the number of differentially expressed genes was decreased after the 35 minute time point at 14 °C and 28 °C (Figure S3, Figure S10), indicating that adaptation started after 35 minutes in P. oligofermentans at 14 °C and 28 °C. Results of the GO enrichment analysis showed that RNA methylation, RNA processing, translation, ribosome biogenesis and cell division GO terms were enriched for upregulated genes at all cold temperatures (Fig. 9). Enrichment of the other modification terms, such as tRNA modification, rRNA processing, and methylation GO terms was similarly detected at cold temperatures implying a role for methylation and modification processes in cold adaptation of P. oligofermentans. Protein folding and de novo UMP biosynthesis were the most enriched GO terms for the upregulated genes at 28 °C (Fig. 9), alongside with the significant enrichment of ATP synthesis and cell division related GO terms for downregulated genes at 28 °C (Figure S11).

Cold-shock, heat-shock and stress related genes in P. oligofermentans

Two cold-shock protein genes were annotated in P. oligofermentans and interestingly gene expression level reactions to cold temperatures of these genes were opposite. Cold-shock protein gene cspA was significantly upregulated at cold temperatures, while cold-shock protein gene cspD was downregulated (Fig. 10b). Several RNase genes (Ribonuclease J1, J2, HII, Z, R, M5, P) and DEAD-box RNA helicase genes (cshA, cshB) were also significantly upregulated at cold temperatures (Fig. 10a,10c). Cold-induced nusA-IF-2 operon was found in P. oligofermentans and the operon was upregulated at cold temperatures, as well as the other transcription initiation factors, IF-1 and IF-3 (Fig. 10d).

Heat-shock related genes, such as heat-inducible transcription repressor hrcA, chaperone genes (groS, groL and dnaK) and Clp protease genes (clpC, clpA), were upregulated at 28 °C at 35 minute time point (Fig. 10e) and downregulated at cold temperatures.

Five universal stress protein (Usp) genes were annotated in P. oligofermentans. Similarly to Le. gelidum and Lc. piscium, downregulation of stress genes was detected at cold temperatures. Only UspA4 was significantly upregulated at 14 °C at 185 minute time point. The UspA4 and PspC family stress protein genes were upregulated at 28 °C (Fig. 10f).

Gene network inference in P. oligofermentans

Based on gene network inference and clustering analysis, 92 clusters contained more than two genes (Table S14). Cold-shock response genes were seen in four different clusters (Cl4, 6, 7, 32), but these clusters were significantly linked. Cl4 (including cshB) was enriched for ribosome biogenesis GO term, Cl6 (including cshA) was enriched for translation GO term, Cl7 (including ribonuclease HII and M5 protein encoding genes) was enriched for rRNA processing GO term, and Cl32 (including ribonuclease R gene) was enriched for regulation of cell shape GO term (Fig. 11). Pseudouridine gene and several number of methylation genes were found in cold-related clusters, which indicated there was a strong interaction between pseudouridine genes, methylation genes, and cold-shock response genes, and suggests methylation and pseudouridine played a role in cold adaptation in P. oligofermentans. Two-component system yycFG gene was also seen in Cl4, which suggests that the two-component system yycFG gene was responsible for cold sensing. We predicted that ribosome related genes and ribosomal protein genes, such as ribosome biogenesis, ribosomal silencing factor RsfS, 50S ribosomal protein L33, 50S ribosomal protein L20, and 30S ribosomal protein S12, were also part of cold-shock response in P. oligofermentans due to strong interaction with cold-shock response genes. Heat-shock related genes and putative TetR family transcriptional regulator gene were clustered within Cl3, indicating the potential role of TetR in heat-shock genes regulation.

de novo transcription binding site prediction in P. oligofermentans

The most enriched known transcription factor binding site motif across all the upstream regions of P. oligofermentans genes was the MalT motif, while rpoD17, CcpA, GalS, GalR, Crp, and SigH motifs were also significantly enriched (Table S15). CspA-binding site motif search across all upstream regions showed that 306 upstream regions contained the 5'-ATTGG-3' motif. This motif was detected in the upstream regions of several cold-induced genes, such as ribonuclease genes, RNA methyltransferase genes, and few heat-induced genes, such as the clpA gene (Table S16).

de novo motif discovery of upregulated genes at different temperatures showed that most of the discovered motifs were similar at each time point. Although for some motifs, tomtom database match analysis showed that they matched with the MalT motif from PRODORIC database [33], these motifs were more likely Shine-Dalgarno sequence motifs of ribosomal binding sites (Table S17). The large number of upstream regions that contributed to the construction of the motif also supported that. Besides the Shine-Dalgarno sequence motif, four discovered motifs had significant scores, but these motifs did not match with the TFBS databases (Table S17).

Statistically significant motif was found for nine clusters based on de novo motif discovery. The discovered motifs for Cl10 and Cl30 were matched with the CcpA motif. In addition, discovered motifs for Cl2 and Cl28 matched with the GalR motif and the GalRS motif, respectively (Table S18). The other discovered motifs were not matched with any motif in the TFBS database.

Pathway enrichment and changes of metabolism at different temperatures in P. oligofermentans

KEGG pathway enrichment analysis showed that peptidoglycan biosynthesis, ribosome and vancomycin resistance KEGG pathway terms were enriched at all cold temperatures for upregulated genes (Fig. 12). This indicated that ribosome and cell wall activities played a role in cold adaptation in P. oligofermentans. Enrichment of aminoacyl-tRNA biosynthesis KEGG term at 0 and 4 °C suggested that production of aminoacyl-tRNA is part of the cold-shock response. Only a few KEGG pathways were enriched for upregulated genes at 28 °C: pyrimidine metabolism was enriched at 5 minute time point and penicillin biosynthesis term was enriched at 185 minute time point (Fig. 12). Carbon and central metabolism terms were enriched for downregulated genes at cold temperatures (Figure S12).

Activity of genes linked with food spoilage

We also specifically checked the gene expression reaction of spoilage genes at both cold- and heat-shock temperatures in the three species. Mostly spoilage-related genes were upregulated at cold temperatures in Le. gelidum (Fig. 13). In addition, slime-related eps genes were upregulated in both Lc. piscium and P. oligofermentans.

Validating RNA-Seq Results With ddPCR

To validate the differential expression results obtained using RNA-seq data, we measured expression levels of selected eight genes (mapA, nagB, pfl, fruK, infC, ftsQ, infB, and 16S rRNA) using droplet digital PCR (ddPCR). Comparison of log₂ fold changes of RNA-seq data and ddPCR data revealed significant correlation (Pearson's correlation coefficient = 0,94) between the two methods (Fig. 14).

Cold-shock and heat-shock responses are both very old systems and share features among bacteria and also with multicellular organisms to some extent. We were interested in the transcription-level reactions of Leuconostoc gelidum, Lactococcus piscium and Paucilactobacillus oligofermentans to both cold and heat. The experiments conducted will help us to understand how these reactions are organized on a level of genomes unraveling LAB function also in the temperatures relevant for food storage.

Bacteria can harbour several CSP gene paralogues, nevertheless not all CSP genes are cold-induced. For example, only four (cspA, cspB, cspG and cspI) of the total nine CSP genes are cold-induced in E.coli [35]. Similarly, it has been shown that four (cspA, cspB, cspC and cspD) of the total five CSP genes are cold-induced in mesophilic lactic acid bacteria Lactococcus lactis [36]. We showed that Le. gelidum and P. oligofermentans harboured two CSP genes (cspA and cspB/cspD), and interestingly Lc. piscium harboured only one CSP gene, cspA. In all three species only cspA was cold-induced, i.e. significant upregulation was observed at cold temperatures based on RNA-sEq. This indicates that cspA was the main responsible CSP gene in these strains for cold adaptation, while cspB/cspD genes may not play a role in cold adaptation in Le. gelidum and P. oligofermentans. DEAD-box RNA helicase genes participate in degradation and unwinding of RNA, and ribosome biogenesis during cold-shock [37]. Our RNA-seq results revealed that DEAD-box RNA helicase genes (cshA in Le. gelidum, cshA and cshB in Lc. piscium and P. oligofermentans) were one of the most significantly upregulated genes in all three species at cold temperatures (Table S1, S7, S13). This suggests that DEAD-box RNA helicase genes have a critical role in these three species at cold temperatures. In addition to DEAD-box RNA helicase, we observed that RNase R and RNase J were the only cold-induced RNase genes in Le. gelidum, while several RNase genes were cold-induced in Lc. piscium and P. oligofermentans (Fig. 6c, Fig. 10c). Though RNase R was reported as the major RNase gene in E. coli [24], upregulation of several RNase genes in Lc. piscium and P. oligofermentans indicates that also other RNase genes took part in cold adaptation.

Translation regulation through translation initiation factors (IF-1, IF-2, IF-3) at cold temperatures has been shown by previous studies [38]. Especially IF-3 significantly favors translation of CSPs [39]. In all three species IF-3 gene was significantly upregulated at cold temperatures indicating similar translation regulation in these species. In addition, IF-2 has been shown to be located in the same operon with nusA in E. coli, and nusA-IF-2 operon affects ribosome maturation at cold temperatures [30]. We observed that in all three species IF-2 was also in the same operon with nusA. The significant upregulation of nusA-IF-2 in Le. gelidum and P. oligofermentans at cold temperatures showed nusA-IF-2 operon can be important for cold adaptation in these species. However, no upregulation of the nusA-IF-2 operon was observed in Lc. piscium.

Besides commonly known cold-shock response genes, we also observed several other cold-induced genes in all three species during cold-shock treatments. The functions of these genes were summarised using GO enrichment and KEGG enrichment analysis (Fig. 1, 4 (Le. gelidum), 5, 8 (Lc. piscium), 9 and 12 (P. oligofermentans)). Ribosome biogenesis, methylation, RNA processing/RNA catabolic process, and fatty acid/lipid metabolism GO and KEGG terms were enriched for upregulated genes at cold temperatures, indicating a similar reaction to cold by the three species. Not all significantly differentially expressed genes need to be related to cold-shock response, since also growth was significantly affected by temperature change. High number of differentially expressed genes in Lc. piscium may be partially explained by different growth phases as was shown in our previous transcriptomics study [14]. To directly identify the temperature change shock response genes, we took advantage of our complex study design (time-series study with several different conditions). It gave us an opportunity to perform not only differential expression gene analysis but also gene network inference analysis and clustering based on gene expression profiles [40]. Clustered genes and linked clusters by gene network are expected to be functionally related [41, 42, 43], and inspecting clusters allowed us to predict genes directly related to the temperature change response.

Gene network inference and clustering results suggests that ribosomal protein and RNA methyltransferase genes were part of cold-shock response in all three species since 50S ribosomal protein L33 gene, RNA/rRNA methyltransferase genes and ABC transporter ATP-binding protein gene were clustered with known cold-shock response genes (Fig. 15).

Specifically, one ABC transporter gene (CBL90879.1, CEN27266.1, CUS25589.1 in Le. gelidum, Lc. piscium and P. oligofermentans, respectively) and one multidrug efflux MFS transporter gene (CBL92085.1, CEN28822.1, CUS25692.1 in Le. gelidum, Lc. piscium and P. oligofermentans, respectively) were clustered with cold-shock response genes and significantly upregulated during cold-shock treatment in all three species, indicating that functionality of these genes was also related to cold-shock response. Both ABC and MFS transporters can play a role in stress resistance, especially in antibiotics resistance [44]. However, previous studies related to psychrophilic organisms also support that ABC and MFS transporters play a role in cold adaptation [45, 46]. Upregulation of multidrug ABC transporter gene during cold temperature was reported for psychrophilic gram-negative bacteria Flavobacterium psychrophilum [45]. Similarly, upregulation of a MFS transporter gene during cold temperature was seen in psychrophilic fungus Mrakia psychrophila [46].

Ribosomal RNA methylation as a response to stress conditions has been shown before. For example, deletion of ksgA (encoding RNA small subunit methyltransferase A) in E. coli causes cold sensitivity and changes in RNA biogenesis [47]. A large number of upregulated RNA methyltransferase genes (Figure S13) at cold temperatures and clustering with cold-shock response genes revealed that RNA methylation plays a role in cold adaptation of the three species studied as a part of post-transcriptional regulation.

Although tRNA methyltransferase and cold-shock genes interaction was only detected in Lc. piscium and P. oligofermentans, tRNA modification gene, tRNA pseudouridine synthase, and tRNA ligase genes clustered with cold-shock genes in Le. gelidum also (Fig. 15). This indicated that tRNA modification was part of cold adaptation in all three species. A tRNAome study showed that tRNA abundance and tRNA modification significantly changed under stress in Lactococcus lactis [48]. It is interesting that tRNA modification is required for cold adaptation in extreme-thermophilic bacteria [49]. In addition, it has been shown that tRNA modification during environmental stress is critical [50]. The gene trmD (encoding guanine-N(1)-methyltransferase) is required for multi-drug resistance in E. coli, since it provides a control mechanism to translation of membrane proteins [50]. It is possible that tRNA modification in the three lactic acid bacteria studied provides a similar control system for cold-shock response genes, promoting their translation.

Two-component sensor histidine kinase activity of DesK-DesR at cold temperature has been shown previously in Bacillus subtilis [51]. The homolog of DesK-DesR was not found in the species studied. However, several other sensor genes were upregulated during cold-shock in all three species (Figure S14). Genes that encode Two-component system WalR (yycF) and sensor kinase WalK (yycG) had similar expression profiles with cold-shock response genes in Le. gelidum and P. oligofermentans (Fig. 15), implying these sensor kinase genes may be linked to cold-shock response in both strains. Previously, sensor kinase WalK has not been reported as a temperature sensor, but it has been shown to have a role in regulating cell wall modifications and peptidoglycan biosynthesis [52]. We suggest that Le. gelidum and P. oligofermentans used WalK to regulate cell wall changes in order to adjust for cold temperature. Interestingly, the homologs of these genes (walK-walR) in Lc. piscium were not upregulated at cold temperatures. In addition, hpk4 gene (CBL92274.1), encoding a two-component sensor histidine kinase was clustered with cold-shock response genes in Le. gelidum, and it was one of the few upregulated genes on early time points at cold temperatures (Figure S14, Table S1). Thus, we believe that hpk4 is a potential cold-sensor gene and could play a role in regulation of cold-shock response genes in Le. gelidum.

In all three species, well known heat-shock response genes, such as chaperone genes, Clp protease genes, and heat-inducible transcriptional repressor hrcA [26], were significantly upregulated at 28 °C (Fig. 2e, 6e, 10e). The main activities of heat-shock response genes are maintaining the proper protein-folding mechanism in cells [26]. Unsurprisingly, protein folding was the only common GO term enriched for upregulated genes in all three species at 28 °C. As well as protein-folding, in Le. gelidum central metabolism GO terms were enriched for upregulated genes at 28 °C, indicating that central metabolism benefited from the temperature increase. We also observed that glycerol and glycerololipid metabolism related genes were upregulated at 28 °C in Lc. piscium. Likely one of the heat-shock responses in Lc. piscium was to change formation of glycerol-based membrane lipids [53]. Notably, a large number of genes clustering with the known cold-shock response genes, potentially signifying the difference in cold-shock adaptation mechanisms among different LAB species. However, heat-shock adaptation mechanisms seem to be highly conservative in LAB and were provided by only the known heat-shock response genes. We also noticed that ctsR gene, a negative regulator of heat-shock Clp-family gene expression [54], was not harboured by Le. gelidum, while it was upregulated at heat-shock temperature in Lc. piscium and P. oligofermentans. Since ctsR is a negative regulator of Clp genes and these genes are required for heat-shock adaptation [54], lack of ctsR gene might be an advantage for Le. gelidum at heat-shock temperatures.

The LAB species studied here are found in cold-stored modified-atmosphere-packaged (MAP) food products and they are the main reason why these food products get spoiled [3, 4]. Therefore, we were interested in the gene expression reactions of spoilage-related genes [12, 14] at different temperatures, especially at cold temperatures. We showed that most of the spoilage-related genes were upregulated at cold temperatures in Le. gelidum. In addition, slime-related genes were upregulated at cold temperatures in Lc. piscium and P. oligofermentans. Thus, we predicted that the level of spoilage activities at cold temperatures were higher for these strains compared to 25 °C.

We were also specifically interested in transcription factor binding site motifs. We showed that CcpA binding site motif was enriched for upregulated genes at both cold- and heat-shock temperatures. It has been shown that CcpA regulon controls carbon metabolism in bacteria [55]. We observed that expression of carbon metabolism genes was significantly affected by temperature change in Le. gelidum. Therefore, enrichment of CcpA motif could be expected for upregulated genes at both cold and warm temperatures. In addition, we discovered motif in the upstream regions of heat-shock genes (from Cl3), which was significantly similar to the HrcA binding site from RegPrecise database [32]. Similarly, analysis of upstream regions for the heat-shock genes clustered in Cl3 in Lc. piscium allowed us to predict binding site motif for heat-shock genes, which was significantly similar to CtsR binding site from RegPrecise database [32]. The motif discovery analysis did not show significant motifs for heat-shock genes in P. oligofermentans. However, we detected motifs for carbon metabolism clusters, which were significantly similar to the CcpA binding site motif. This shows ccpA played a role in metabolism response to temperature change, similar to Le. gelidum.

Here, we profiled gene expression changes of three psychrotrophic lactic acid bacteria at three time points and five different temperatures. The results provided an important understanding of the cold- and heat-shock adaptation mechanisms in psychrotrophic LAB. Besides known cold-shock response genes, we highlighted the importance of rRNA and tRNA modification during cold-shock response in psychrotrophic LAB as a part of post-transcriptional regulation. In addition, we showed that ABC and efflux MFS transporter genes were involved in cold-shock response. These results might give new directions to cold-shock studies in psychrotrophic organisms. Transcriptomic profiles during heat-shock treatment revealed that protein folding was the dominant reaction for psychrotrophic LAB. Upregulation of spoilage-related genes at cold temperatures suggests the level of spoilage activities was higher at cold temperatures, especially for Le. gelidum. Further gene inference analysis revealed several candidate genes involved in cold- and heat-shock adaptation and sensing. The candidate genes identified using gene inference analysis will be beneficial for further studies. Finally, we showed binding site motifs for heat-shock genes in Le. gelidum and Lc. piscium using de novo motif search analysis.

Experimental Setup and Sampling

Three psychrophilic lactic acid bacteria (LAB) species: Leuconostoc gelidum subsp. gasicomitatum LMG 18811, Lactococcus piscium MKFS47, Paucilactobacillus oligofermentans DSM 15707, were used for the experiment. Frozen glycerol stock solution (1.8 ml) of a bacterium was inoculated to 9 ml of de Man-Rogosa-Sharpe medium without acetate (MRS-Ac) and incubated overnight at 25 °C. One µl of the culture was streaked on a MRS-Ac agar plate and grown in anaerobic conditions over three nights at 25 °C. One colony was then transferred to 9 ml of MRS-Ac medium. After 17 h (for Le. gelidum and Lc. piscium) or 31 h (for P. oligofermentans), a 100-µl aliquot was again inoculated to 9 ml of MRS-Ac medium and grown for further 12 h (for Le. gelidum and Lc. piscium) or 24 h (for P. oligofermentans). An aliquot of 90 µl (Le. gelidum and Lc. piscium) or 80 µl (P. oligofermentans) of these precultures was inoculated into 100 ml MRS-Ac and grown over one night (Le. gelidum and Lc. piscium) or over two nights (P. oligofermentans) at 25 °C until the OD₆₀₀ value was about 0.4. The cultures were then diluted 1:10 to have an OD₆₀₀ value 0.04 in a total volume of 100 ml. The obtained cultures were grown at 25°C in four replicates. The growth of these cultures was followed by making OD₆₀₀measurements every hour for 10 h (Le. gelidum and Lc. piscium) or 13 h (P. oligofermentans). After 4h (Lc. piscium), 5h (Le. gelidum), and 9h (P. oligofermentans) of growth, 10-ml aliquots of each culture were transferred to five different temperatures: 0 °C, 4 °C, 14 °C, 25 °C, and 28 °C. The temperatures were selected to include cold-shock temperatures (0 °C, 4 °C), mild cold-shock temperature (14 °C) and a heat-shock temperature (28 °C), while 25 °C was used as the control temperature. For each bacterium, samples for RNA extraction were collected from the 100-ml bottles before aliquoting (time point 0 min), and from the aliquots in different temperatures after 5 min, 35 min and 185 min of growth (Figure S1). Sample volumes were determined using the formula: V = (1/1.5*OD₆₀₀). The same sample volume was used at 0 min, 5 min and 35 min time points, while before sampling at 185 min, OD₆₀₀ measurement was made to adjust the sample volume. Samples were immediately mixed with a 1/10 of a volume of cold ethanol-phenol mixture (10:1) to stop RNase activity. Cells were pelleted by centrifugation at 5000 x g for 3 min, frozen in liquid nitrogen, and stored in -80°C. Altogether 192 samples were collected during the growth experiment

RNA extraction and transcriptome sequencing

RNA extraction was performed with Nucleospin RNA kit (Macherey-Nagel, Düren, Germany) according to manufacturer’s instructions with some modifications in the cell lysis step. Cells were resuspended in 700 µl of buffer RA1 supplemented with 7 µl of β-mercaptoethanol, and then mechanically disrupted using Lysing Matrix E or B tubes (MP Biomedicals, Irvine, CA, USA) and FastPrep 24 tissue homogenizer (MP Biomedicals) at 5.5 m/s for 40 sec. Quality and quantity of RNA extractions were analyzed using the Agilent 2100 Bioanalyzer and RNA 6000 Nano kit (Agilent, Santa Clara, CA, USA).

Ribosomal RNA removal was performed with Ribo-Zero rRNA removal reagent for bacteria (Illumina, San Diego, CA, USA) according to manufacturer’s instructions using 1000 ng of total RNA and a 1/3 volume of kit solutions. rRNA-depleted RNA was purified with RNA Clean&Concentrator -5 kit (Zymo Research, Irvine, California, USA) according to manufacturer’s instructions and eluted in 11 µl of RNA-free water.

QIAseq™ Stranded Total RNA Lib kit (Qiagen, Hilden, Germany) was used for RNA-seq library preparation according to manufacturer’s instructions with some modifications. One third of a volume of kit solutions were used together with 9.7 µl rRNA-depleted RNA. In strand-specific ligation step, 0.3 µM truncated Truseq adapter was used instead of the adapter provided in the kit. Libraries were amplified using half (10 µl) of the purified ligation product in 1 x HF buffer with 0.2 mM dNTPs, 0.6 µM selected [56] dual-index primer, and one unit of Phusion Hot Start II High-Fidelity DNA polymerase (Thermo Fisher Scientific, Waltham, Massachusetts, United States) in a total volume of 50 µl. PCR protocol used was 98°C, 30 sec; 18 x (98°C, 10 sec; 65°C, 30 sec; 72°C, 10 sec); 72, 5 min. Concentrations of amplified libraries were measured with Qubit fluorometer and dsDNA HS assay kit (Invitrogen, Waltham, MA, USA), and size distributions visualized with Fragment Analyzer and High Sensitivity NGS Fragment Analysis kit (Advanced Analytical, Parkersburg, WV, USA). After pooling, libraries were concentrated using Amicon Ultra 100K columns (Millipore, Burlington, MA, USA), and purified twice with 0.9 x AMPure XP beads (Beckman Coulter, Brea, CA, USA). For Le. gelidum and P. oligofermentans libraries, size selection of 300 - 600 bp fragments was performed using BluePippin and 2% agarose gel cassette (Sage Science, Beverly, MA, USA). With Lc. piscium libraries PEG (8% - 8.5%) size selection was used. Illumina NextSeq 500 was used to sequence the resulting RNA-seq libraries single-end 75 bp.

Sequencing data processing

Adaptor trimming and quality filtering of the reads was done using Trimmomatic v0.36 [57] using these parameters “TruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:3:20 HEADCROP:10 MINLEN:30”. SortMeRNA v2.1 [58] was used for ribosomal RNA filtering. Filtered reads were mapped to genomes using Bowtie2 v2.3.4.3 [59] and were sorted using Samtools v1.8 [60]. Read counts for genes were produced using HTSeq v0.11.0 [61] with union mode. DESeq2 v1.22.2 R package [62] was used for differential gene expression analysis. The significance threshold for differential gene expression was set for P-adjusted <= 0.05 and |log2FoldChange| => 1. PANNZER2 [63] was used for the functional and GO term annotation. The enrichment analysis for GO terms was done using topGO v2.36.0 R package [64], the enrichment threshold was set for p-value <= 0.05. KEGG pathway enrichment was done using enrichKEGG function of clusterProfiler v3.14.3 R package [65]. KEGG annotation of Le. gelidum and Lc. piscium was obtained from the KEGG database [66]. Since P. oligofermentans was not found in the KEGG database, KAAS (KEGG Automatic Annotation Server) webserver [67] was used to create KEGG annotation and KEGG pathways of P. oligofermentans. We also used Pathway Tools v19.0 [68] for pathway construction and enrichment.

Gene Network inference analysis and Motif Discovery

Gene network inference was done using seidr toolkit [40] using 11 gene network inference methods; CLR, GENIE3, Aracne2, Pearson, Spearman, NARROMI, TIGRESS, PCor, PLSNET, llr-ensemble, and el-ensemble. Their results were aggregated into a network using seidr toolkit. Hard threshold was chosen following the recommendations from the developer. The network clustered using Infomap [69] and clusters visualized using the Map generator applet from MapEquation [70].

For the motif based sequence analysis MEME suite v5.0.5 [71] was used. The upstream regions had minimal length of 50 and up to 300 nucleotides were extracted using python script (https://github.com/peterthorpe5/intergenic_regions). Motif discovery of upstream regions was done using MEME, discovered motifs were search against transcription factor binding site databases, such as CollecTF [72], PRODORIC [33], RegTransBase [73], RegPrecise [32], DPINTERACT [34], and Swiss Regulon [74], using Tomtom [75]. The listed databases were downloaded from MEME suite as meme database format, except RegPrecise. For RegPrecise, all transcription factor binding site motifs for Lactobacillaceae were downloaded and converted to meme database using sites2meme script from MEME suite v5.0.5 [71]. Ame [76] was used to test whether motifs in transcription binding site databases enriched in upstream regions of our reference genomes.

ddPCR validation

Three out of four replicate samples grown at 4 °C, 14 °C, 25 °C, and 28 °C, and collected at 185 min were analyzed by ddPCR to verify the RNA-seq results. Primers and the protocol used were previously published by Andreevskaya et al. [16]. For Le. gelidum, expression of mapA, nagB, ftsQ, infB and 16S rRNA genes, for Lc. piscium, expression pfl, fruK, ftsQ, infB, and 16S rRNA genes , and for P. oligofermentans, expression of ccpA, surf, infC, ftsQ, infB, and 16S rRNA genes were analyzed. For each species, the reference genes used for normalization of gene copy numbers were the ones that showed the most stable expression levels (according to RNA-seq) among the samples studied: infB, ftsQ, and 16S rRNA gene (Le. gelidum and Lc. piscium), or ccpA, surf, and 16S rRNA gene (P. oligofermentans).

cold-shock protein (CSP), lactic acid bacteria (LAB), ribonuclease (RNase), GO (Gene Ontology), Cluster (Cl), transcription factor binding site (TFBS), KEGG (Kyoto Encyclopedia of Genes and Genomes), droplet digital PCR (ddPCR), initiation factor (IF)

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

All sequencing data have been deposited in the European Nucleotide Archive (ENA) under accession code PRJEB38386.

Competing interests

The authors declare that they have no competing interests.

Funding

The study was funded by the Academy of Finland project 307856 and 323576 to PA. ICD supported by the MBDP doctoral program. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Authors' contributions

PA conceived and designed the study. AY and AAP collected the samples, performed the temperature shock treatment, RNA extraction and sequencing preparation. LP organized NGS assays. ICD performed the bioinformatics analyses and drafted the manuscript with input from AY, PA, and SB. All authors have read, commented, and approved the final manuscript.

Acknowledgements

We want to thank Eeva-Marja Turkki and Päivi Laamanen for performing the NGS procedures for this project. We wish to acknowledge CSC – IT Center for Science, Finland, for computational resources.

Carr FJ, Chill D, Maida N. The Lactic Acid Bacteria: A literature survey. Crit Rev Microbiol. 2002;28:281–370.
Singh VP. Recent approaches in food bio-preservation - a review. Open Vet J. 2018;8:104–11.
Remenant B, Jaffrès E, Dousset X, Pilet M-F, Zagorec M. Bacterial spoilers of food: Behavior, fitness and functional properties. Food Microbiol. 2015;45:45–53.
Nieminen TT, Koskinen K, Laine P, Hultman J, Säde E, Paulin L, et al. Comparison of microbial communities in marinated and unmarinated broiler meat by metagenomics. Int J Food Microbiol. 2012;157:142–9.
Ercolini D, Ferrocino I, Nasi A, Ndagijimana M, Vernocchi P, La Storia A, et al. Monitoring of Microbial Metabolites and Bacterial Diversity in Beef Stored under Different Packaging Conditions. Appl Environ Microbiol. 2011;77:7372–81.
Björkroth KJ, Geisen R, Schillinger U, Weiss N, Vos PD, Holzapfel WH, et al. Characterization of Leuconostoc gasicomitatum sp. nov., associated with spoiled raw tomato-marinated broiler meat strips packaged under modified-atmosphere conditions. Appl Environ Microbiol. 2000;66:3764–72.
Vihavainen EJ, Björkroth KJ. Spoilage of value-added, high-oxygen modified-atmosphere packaged raw beef steaks by Leuconostoc gasicomitatum and Leuconostoc gelidum. Int J Food Microbiol. 2007;119:340–5.
Rahkila R, Nieminen T, Johansson P, Säde E, Björkroth J. Characterization and evaluation of the spoilage potential of Lactococcus piscium isolates from modified atmosphere packaged meat. Int J Food Microbiol. 2012;156:50–9.
Jääskeläinen E, Johansson P, Kostiainen O, Nieminen T, Schmidt G, Somervuo P, et al. Significance of heme-based respiration in meat spoilage caused by Leuconostoc gasicomitatum. Appl Environ Microbiol. 2013;79:1078–85.
Doulgeraki AI, Ercolini D, Villani F, Nychas G-JE. Spoilage microbiota associated to the storage of raw meat in different conditions. Int J Food Microbiol. 2012;157:130–41.
Koort J, Murros A, Coenye T, Eerola S, Vandamme P, Sukura A, et al. Lactobacillus oligofermentans sp. nov., associated with spoilage of modified-atmosphere-packaged poultry products. Appl Environ Microbiol. 2005;71:4400–6.
Johansson P, Paulin L, Säde E, Salovuori N, Alatalo ER, Björkroth KJ, et al. Genome sequence of a food spoilage lactic acid bacterium, Leuconostoc gasicomitatum LMG 18811T, in association with specific spoilage reactions. Appl Environ Microbiol. 2011;77:4344–51.
Andreevskaya M, Hultman J, Johansson P, Laine P, Paulin L, Auvinen P, et al. Complete genome sequence of Leuconostoc gelidum subsp. gasicomitatum KG16-1, isolated from vacuum-packaged vegetable sausages. Standards in Genomic Sciences. 2016;11:40.
Andreevskaya M, Johansson P, Laine P, Smolander O-P, Sonck M, Rahkila R, et al. Genome sequence and transcriptome analysis of meat-spoilage-associated lactic acid bacterium Lactococcus piscium MKFS47. Appl Environ Microbiol. 2015;81:3800–11.
Andreevskaya M, Johansson P, Jääskeläinen E, Rämö T, Ritari J, Paulin L, et al. Lactobacillus oligofermentans glucose, ribose and xylose transcriptomes show higher similarity between glucose and xylose catabolism-induced responses in the early exponential growth phase. BMC Genom. 2016;17:539.
Andreevskaya M, Jääskeläinen E, Johansson P, Ylinen A, Paulin L, Björkroth J, et al. Food spoilage-associated Leuconostoc, Lactococcus, and Lactobacillus species display different survival strategies in response to competition. Appl Environ Microbiol. 2018;84. doi:10.1128/AEM.00554-18.
Zeikus JG. Thermophilic bacteria: Ecology, physiology and technology. Enzyme Microbial Technology. 1979;1:243–52.
D’Amico S, Collins T, Marx J-C, Feller G, Gerday C, Gerday C. Psychrophilic microorganisms: challenges for life. EMBO Rep. 2006;7:385–9.
Morita RY. Psychrophilic bacteria. Microbiol Mol Biol Rev. 1975;39:144–67.
Hébraud M, Potier P. Cold shock response and low temperature adaptation in psychrotrophic bacteria. J Mol Microbiol Biotechnol. 1999;1:211–9.
De Maayer P, Anderson D, Cary C, Cowan DA. Some like it cold: understanding the survival strategies of psychrophiles. EMBO Rep. 2014;15:508–17.
Barria C, Malecki M, Arraiano CM. Bacterial adaptation to cold. Microbiology. 2013;159 Pt_12:2437–43.
Tribelli PM, López NI. Reporting key features in cold-adapted bacteria. Life. 2018;8:8.
Zhang Y, Burkhardt DH, Rouskin S, Li G-W, Weissman JS, Gross CA. A Stress response that monitors and regulates mRNA structure is central to cold shock adaptation. Mol Cell. 2018;70:274–86.e7.
Phadtare S. Recent developments in bacterial cold-shock response. Curr Issues Mol Biol. 2004;6:125–36.
Mogk A, Tomoyasu T, Goloubinoff P, Rüdiger S, Röder D, Langen H, et al. Identification of thermolabile Escherichia coli proteins: prevention and reversion of aggregation by DnaK and ClpB. The EMBO Journal. 1999;18:6934–49.
Wouters JA, Rombouts FM, Kuipers OP, de Vos WM, Abee T. The role of cold-shock proteins in low-temperature adaptation of food-related bacteria. Syst Appl Microbiol. 2000;23:165–73.
Kim WS, Dunn NW. Identification of a cold shock gene in lactic acid bacteria and the effect of cold shock on Cryotolerance. Curr Microbiol. 1997;35:59–63.
Varmanen P, Savijoki K. Responses of lactic Acid Bacteria to heat stress. In: Tsakalidou E, Papadimitriou K, editors. Stress responses of lactic acid bacteria. Boston: Springer US; 2011. pp. 55–66. doi:10.1007/978-0-387-92771-8_3.
Brandi A, Giangrossi M, Paoloni S, Spurio R, Giuliodori AM, Pon CL, et al. Transcriptional and post-transcriptional events trigger de novo infB expression in cold stressed Escherichia coli. Nucleic Acids Res. 2019;47:4638–51.
Graumann P, Marahiel MA. The major cold shock protein of Bacillus subtilis CspB binds with high affinity to the ATTGG- and CCAAT sequences in single stranded oligonucleotides. FEBS Lett. 1994;338:157–60.
Novichkov PS, Kazakov AE, Ravcheev DA, Leyn SA, Kovaleva GY, Sutormin RA, et al. RegPrecise 3.0 – A resource for genome-scale exploration of transcriptional regulation in bacteria. BMC Genom. 2013;14:745.
Münch R, Hiller K, Barg H, Heldt D, Linz S, Wingender E, et al. PRODORIC: prokaryotic database of gene regulation. Nucleic Acids Res. 2003;31:266–9.
Robison K, McGuire AM, Church GM. A comprehensive library of DNA-binding site matrices for 55 proteins applied to the complete Escherichia coli K-12 genome. J Mol Biol. 1998;284:241–54.
Yamanaka K, Fang L, Inouye M. The CspA family in Escherichia coli : multiple gene duplication for stress adaptation. Mol Microbiol. 1998;27:247–55.
Woufers JA, sander J-W, Kok J, de Vos WM, Kuipers OP, Abee T. Clustered organization and transcriptional analysis of a family of five csp genes of Lactococcus lactis MGl363. Microbiology. 1998;144:2885–93.
Hunger K, Beckering CL, Wiegeshoff F, Graumann PL, Marahiel MA. Cold-Induced putative DEAD Box RNA helicases CshA and CshB are essential for cold adaptation and interact with cold shock protein B in Bacillus subtilis. J Bacteriol. 2006;188:240–8.
Giuliodori AM, Brandi A, Gualerzi CO, Pon CL. Preferential translation of cold-shock mRNAs during cold adaptation. RNA. 2004;10:265–76.
Giuliodori AM, Brandi A, Giangrossi M, Gualerzi CO, Pon CL. Cold-stress-induced de novo expression of infC and role of IF3 in cold-shock translational bias. RNA. 2007;13:1355–65.
Schiffthaler B, Serrano A, Street N, Delhomme N. Seidr: a gene meta-network calculation toolkit. bioRxiv. 2019;:250696.
Luo F, Yang Y, Zhong J, Gao H, Khan L, Thompson DK, et al. Constructing gene co-expression networks and predicting functions of unknown genes by random matrix theory. BMC Bioinformatics. 2007;8:299.
van Dam S, Võsa U, van der Graaf A, Franke L, de Magalhães JP. Gene co-expression analysis for functional classification and gene–disease predictions. Brief Bioinform. 2018;19:575–92.
Schlitt T, Palin K, Rung J, Dietmann S, Lappe M, Ukkonen E, et al. From gene networks to gene function. Genome Res. 2003;13:2568–76.
Greene NP, Kaplan E, Crow A, Koronakis V. Antibiotic resistance mediated by the MacB ABC transporter family: A structural and functional perspective. Front Microbiol. 2018;9. doi:10.3389/fmicb.2018.00950.
Hesami S, Metcalf DS, Lumsden JS, MacInnes JI. Identification of cold-temperature-regulated genes in Flavobacterium psychrophilum. Appl Environ Microbiol. 2011;77:1593–600.
Su Y, Jiang X, Wu W, Wang M, Hamid MI, Xiang M, et al. Genomic, transcriptomic, and proteomic analysis provide insights into the cold adaptation mechanism of the obligate psychrophilic fungus Mrakia psychrophila. G3 (Bethesda). 2016;6:3603–13.
Connolly K, Rife JP, Culver G. Mechanistic insight into the ribosome biogenesis functions of the ancient protein KsgA. Mol Microbiol. 2008;70:1062–75.
Puri P, Wetzel C, Saffert P, W.Gaston K, Russell SP, Varela JAC, et al. Systematic identification of tRNAome and its dynamics in Lactococcus lactis. Mol Microbiol. 2014;93:944–56.
Ishida K, Kunibayashi T, Tomikawa C, Ochi A, Kanai T, Hirata A, et al. Pseudouridine at position 55 in tRNA controls the contents of other modified nucleotides for low-temperature adaptation in the extreme-thermophilic eubacterium Thermus thermophilus. Nucleic Acids Res. 2011;39:2304–18.
Masuda I, Matsubara R, Christian T, Rojas ER, Yadavalli SS, Zhang L, et al. tRNA methylation is a global determinant of bacterial multi-drug resistance. Cell Systems. 2019;8:302–14.e8.
Aguilar PS, Hernandez-Arriaga AM, Cybulski LE, Erazo AC, de Mendoza D. Molecular basis of thermosensing: a two-component signal transduction thermometer in Bacillus subtilis. EMBO J. 2001;20:1681–91.
Dubrac S, Boneca IG, Poupel O, Msadek T. New insights into the WalK/WalR (YycG/YycF) essential signal transduction pathway reveal a major role in controlling cell wall metabolism and biofilm formation in Staphylococcus aureus. J Bacteriol. 2007;189:8257–69.
Geiger O, Sohlenkamp C, López-Lara IM. Formation of bacterial glycerol-based membrane lipids: Pathways, enzymes, and reactions. In: Geiger O, editor. Biogenesis of fatty acids, lipids and membranes. Cham: Springer International Publishing; 2019. pp. 87–107. doi:10.1007/978-3-319-50430-8_8.
Varmanen P, Ingmer H, Vogensen FK. ctsR of Lactococcus lactis encodes a negative regulator of clp gene expression. Microbiology. 2000;146:1447–55.
Warner JB, Lolkema JS. CcpA-dependent carbon catabolite repression in bacteria. Microbiol Mol Biol Rev. 2003;67:475–90.
Somervuo P, Koskinen P, Mei P, Holm L, Auvinen P, Paulin L. BARCOSEL: a tool for selecting an optimal barcode set for high-throughput sequencing. BMC Bioinformatics. 2018;19:257.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211–7.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Törönen P, Medlar A, Holm L. PANNZER2: a rapid functional annotation web server. Nucleic Acids Res. 2018;46:W84–8.
Alexa A, Rahnenfuhrer J. topGO. topGO: Enrichment analysis for Gene Ontology. Available at: https://bioconductor.org/packages/release/bioc/html/topGO.html.
Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44:D457–62.
Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:W182–5.
Karp PD, Latendresse M, Paley SM, Krummenacker M, Ong QD, Billington R, et al. Pathway Tools version 19.0 update: software for pathway/genome informatics and systems biology. Brief Bioinform. 2016;17:877–90.
Rosvall M, Bergstrom CT. Maps of random walks on complex networks reveal community structure. PNAS. 2008;105:1118–23.
Rosvall M, Axelsson D, Bergstrom CT. The map equation. Eur Phys J Spec Top. 2009;178:13–23.
Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME Suite: tools for motif discovery and searching. Nucleic Acids Res. 2009;37:W202–8.
Kiliç S, White ER, Sagitova DM, Cornish JP, Erill I. CollecTF: a database of experimentally validated transcription factor-binding sites in Bacteria. Nucleic Acids Res. 2014;42:156–60. Database issue:D.
Cipriano MJ, Novichkov PN, Kazakov AE, Rodionov DA, Arkin AP, Gelfand MS, et al. RegTransBase – a database of regulatory sequences and interactions based on literature: a resource for investigating transcriptional regulation in prokaryotes. BMC Genom. 2013;14:213.
Pachkov M, Erb I, Molina N, van Nimwegen E. SwissRegulon: a database of genome-wide annotations of regulatory sites. Nucleic Acids Res. 2007;35 suppl_1:D127–31.
Gupta S, Stamatoyannopoulos JA, Bailey TL, Noble WS. Quantifying similarity between motifs. Genome Biol. 2007;8:R24.
McLeay RC, Bailey TL. Motif Enrichment Analysis: a unified framework and an evaluation on ChIP data. BMC Bioinformatics. 2010;11:165.

Additional File1: Figure S1. Experimental setup and sampling summary. Each of the three strains were grown at 25 °C as four replicates. After collecting the first sample (control; time point 0 min) aliquots were taken for five temperatures; 0 °C, 4 °C, 14 °C, 25 °C, and 28 °C. For each aliquot, samples were collected after 5 minutes, 35 minutes, and 185 minutes.

Additional File2: Figure S2. Growth curve of all three species based on the optical density (OD600) values. The black colored points and line represents growth at 25 °C liquid broth. Sampling times for aliquoting to different temperatures are marked in the figure and. Different colored points represent samples at different temperatures at 185 min; yellow: 0 °C, green: 4°C, blue: 14°C, red: 25°C, and pink: 28°C. Statistically significant (Student’s t-test P-value < 0.05) difference in growth compared to 25°C control aliquot is indicated with *.

Additional File3: Figure S3. Number of differentially expressed genes of all three species at different temperatures. Figure shows the number of differentially expressed genes of three species at 0 °C, 4°C, 14 °C and 28 °C. In general, the numbers of differentially expressed genes were low In the first time point, but increased in the later time points.

Additional File4: Figure S4. Heatmap of log₂ fold changes for all differentially expressed genes in Le. gelidum. Upregulated genes are colored red and downregulated genes are colored blue. Columns represent comparison between different temperatures and control samples.

Additional File5: Table S1. The list of differentially expressed genes in Le. gelidum

Additional File6: Figure S5. Heatmap of enriched GO terms of downregulated genes in Le. gelidum. Enriched GO terms of downregulated genes compared at different temperatures and time points. Comparisons were made against data from the 25 °C control. Red gradient represents enrichment p-value, for which scale is shown at the right corner.

Additional File7: Table S2. The list of gene clusters based on gene network inference analysis in Le. gelidum.

Additional File8: Table S3. The list of known motifs enriched within upstream regions of all genes in Le. gelidum based on analysis of motif enrichment.

Additional File9: Table S4. The list of Le. gelidum genes, which have CspA binding site motif on their upstream regions.

Additional File10: Table S5. The list of motifs with significant p-value discovered based on upstream regions of upregulated genes in Le. gelidum.

Additional File11: Table S6. The list of motifs with significant p-value discovered based on upstream regions of clustered genes in Le. gelidum.

Additional File12: Figure S6. KEGG pathway enrichment for downregulated genes of Le. gelidum. Figure shows heatmap of enriched KEGG pathways for downregulated genes in different temperatures. Enriched KEGG pathways are marked with red. Red gradient represents enrichment p-value, for which scale is shown at the right corner.

Additional File13: Table S7. The list of differentially expressed genes in Lc. piscium

Additional File14: Figure S7. Heatmap of log₂ fold changes for all differentially expressed genes in Lc. piscium. Upregulated genes are colored red and downregulated genes are colored blue. Columns represent comparison between different temperatures and control samples.

Additional File15: Figure S8. Heatmap of enriched GO terms of downregulated genes in Lc. piscium. Enriched GO terms of downregulated genes compared at different temperatures and time points. Comparisons were made against data from the 25 °C control. Red gradient represents enrichment p-value, for which scale is shown at the right corner.

Additional File16: Table S8. The list of gene clusters based on gene network inference analysis in Lc. piscium.

Additional File17: Table S9. The list of known motifs enriched within upstream regions of all genes in Lc. piscium based on analysis of motif enrichment.

Additional File18: Table S10. The list of Lc. piscium genes, which have CspA binding site motif on their upstream regions.

Additional File19: Table S11. The list of motifs with significant p-value discovered based on upstream regions of upregulated genes in Lc. piscium.

Additional File20: Table S12. The list of motifs with significant p-value discovered based on upstream regions of clustered genes in Lc. piscium.

Additional File21: Figure S9. KEGG pathway enrichment for downregulated genes of Lc. piscium. Figure shows heatmap of enriched KEGG pathways for downregulated genes in different temperatures. Enriched KEGG pathways were marked as red. Red gradient represents enrichment p-value, for which scale is shown at the right corner.

Additional File22: Table S13. The list of differentially expressed genes in P. oligofermentans

Additional File23: Figure S10. Heatmap of log₂ fold changes for all differentially expressed genes in P. oligofermentans. Upregulated genes are colored red and downregulated genes are colored blue. Columns represent comparison between different temperatures and control samples.

Additional File24: Figure S11. Heatmap of enriched GO terms of downregulated genes in P. oligofermentans. Enriched GO terms of downregulated genes compared at different temperatures and time points. Comparisons were made against data from the 25 °C control. Red gradient represents enrichment p-value, for which scale is shown at the right corner.

Additional File25: Table S14. The list of gene clusters based on gene network inference analysis in P. oligofermentans.

Additional File26: Table S15. The list of known motifs enriched within upstream regions of all genes in P. oligofermentans based on analysis of motif enrichment.

Additional File27: Table S16. The list of P. oligofermentans genes, which have CspA binding site motif on their upstream regions.

Additional File28: Table S17. The list of motifs with significant p-value discovered based on upstream regions of upregulated genes in P. oligofermentans.

Additional File29: Table S18. The list of motifs with significant p-value discovered based on upstream regions of clustered genes in P. oligofermentans.

Additional File30: Figure S12. KEGG pathway enrichment for downregulated genes of P. oligofermentans. Figure shows heatmap of enriched KEGG pathways for downregulated genes in different temperatures. Enriched KEGG pathways are marked with red. Red gradient represents enrichment p-value, for which scale is shown at the right corner.

Additional File31: Table S19. List of spoilage related genes and log₂ fold changes in the three species.

Additional File32: Figure S13. The log₂ fold change heatmap of methylation-related genes in the three species. The log₂ fold change scale at the right corner.

Additional File33: Figure S14. The log₂ fold change heatmap of sensing/signal-related genes in all three species. The log₂ fold change scale at the right corner.

Download PDF

Journal Publication

published 07 Jan, 2021

Read the published version in BMC Genomics →

Editor assigned by journal
08 Jun, 2020
Submission checks completed at journal
07 Jun, 2020
Editor invited by journal
07 Jun, 2020
First submitted to journal
04 Jun, 2020

You are reading this latest preprint version

Transcriptomic time-series analysis of cold- and heat-shock response in psychrotrophic lactic acid bacteria

Status:

Journal Publication

Version 1

Abstract

Background:

Results:

Conclusions:

Figures

Background

Results

Bacterial growth

Activity of genes linked with food spoilage

Validating RNA-Seq Results With ddPCR

Discussion

Conclusions

Methods

Experimental Setup and Sampling

RNA extraction and transcriptome sequencing

Sequencing data processing

Gene Network inference analysis and Motif Discovery

ddPCR validation

List Of Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and material

Competing interests

Funding

Authors' contributions

Acknowledgements

References

Additional Files

Supplementary Files

Status:

Journal Publication

Version 1