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 log2 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  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 , 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 . 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. Log2 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  and at 28 °C 185 minute time point, discovered motif matched with rpoD17 motif from DPInteract database  (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 , 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).