Metabolic robustness to growth temperature of a cold adapted bacterium

Microbial communities experience continuous environmental changes, among which temperature fluctuations are arguably the most impacting. This is particularly important considering the ongoing global warming but also in the “simpler” context of seasonal variability of sea-surface temperature. Understanding how microorganisms react at the cellular level can improve our understanding of possible adaptations of microbial communities to a changing environment. In this work, we investigated which are the mechanisms through which metabolic homeostasis is maintained in a cold-adapted bacterium during growth at temperatures that differ widely (15 and 0°C). We have quantified its intracellular and extracellular central metabolomes together with changes occurring at the transcriptomic level in the same growth conditions. This information was then used to contextualize a genome-scale metabolic reconstruction and to provide a systemic understanding of cellular adaptation to growth at two different temperatures. Our findings indicate a strong metabolic robustness at the level of the main central metabolites, counteracted by a relatively deep transcriptomic reprogramming that includes changes in gene expression of hundreds of metabolic genes. We interpret this as a transcriptomic buffering of cellular metabolism, able to produce overlapping metabolic phenotypes despite the wide temperature gap. Moreover, we show that metabolic adaptation seems to be mostly played at the level of few key intermediates (e.g. phosphoenolpyruvate) and in the cross-talk between the main central metabolic pathways. Overall, our findings reveal a complex interplay at gene expression level that contributes to the robustness/resilience of core metabolism, also promoting the leveraging of state-of-the-art multi-disciplinary approaches to fully comprehend molecular adaptations to environmental fluctuations.


Introduction
Microorganisms are able to colonize virtually every environmental niche on Earth (Rothschild and Mancinelli, 2001). They have adapted for millions of years prospering under conditions such extreme as water boiling or freezing points, high radiation, acidic or alkaline pH values, heavy metal pollution, and high salinity (Shu and Huang, 2022). Growth temperature, in particular, is one of the environmental parameters that mostly impact the physiology of microorganisms and that is thought to have played a key role in their adaptation, selection and diversification (Hegedűs and Csonka, 2010). Given the geological history of our planet, it is reasonable to think that adaptation to changing temperatures has independently occurred many times in evolution; consequently, there exists a vast array of molecular strategies for this purpose disseminated in the microbial kingdom (Bennett and Lenski, 1997;Saarinen et al., 2018;Walworth et al., 2020). Their characterization is key in this phase of Earth's life, as global change is imposing rapid/drastic modifications in basic environmental parameters (including temperature) that, in turn, will solicitate the activation of such temperature-adaptation related pathways in microbial communities. Understanding which genes get activated or which compounds get secreted in the environment following a future increase of water temperature will help us model and predict the scenarios of microbial communities in a changing environment.
Without necessarily invoking global change, the need to rewire cellular networks in response to temperature shifts is likely a common feature in natural microbial communities. The temperature in the Southern Ocean, for example, is anywhere from -2 to 10°C, since Antarctic water temperature fluctuation responds to the seasonal advance and retreat of sea ice (Auger et al., 2021). Therefore, marine microorganisms are exposed to seasonal oscillations in temperature. Upper-ocean microbes can experience a higher variability in sea-surface temperature and their working temperatures exceed the in situ Eulerian temperature range by up to 10 °C (Doblin and van Sebille, 2016). In general, drifting microbial communities experience large temperature variations and have broader thermal tolerance compared to communities from stable environments (Gilchrist, 1995). Delille at al. suggested that in Antarctic seawater, chemoautotrophic bacterial biomass abundance seems effectively regulated by temperature, whereas its influence on the biomass of heterotrophic bacterial communities is limited (Delille, 1993(Delille, , 1996. Their findings suggest that Antarctic marine heterotrophic bacteria appeared well adapted to temperature fluctuation and are in accordance with physiological studies on isolated strains that clearly show that a large majority of Antarctic heterotrophic marine bacteria have good growth rates for a wide range of temperatures 0 to 20°C (Nedwell and Rutter, 1994). From a cellular viewpoint, previous studies of cold-and heat-adapted microbes have revealed a variety of molecular adaptations that allow their activity and survival under extreme conditions. These initially involve the change in the expression of specific gene sets but, ultimately, the regulatory changes imposed by temperature increase/decrease are mostly implemented at the metabolic level, since this represents, in the words of Prof. Oliver Fiehn, "the ultimate response of biological systems to genetic or environmental changes" (Fiehn, 2002). These variegate adaptation strategies observed in temperature-stressed microorganisms indicate that i) each microorganism may follow a peculiar route to maintain cellular homeostasis when facing temperature fluctuations (Lu et al., 2016;Czajka et al., 2018;Kurdrid et al., 2020) and ii) there exists intense crosstalk between regulatory and metabolic networks in the response thereof (Jozefczuk et al., 2010).
Building upon previous work (Fondi et al., 2015;Mocali et al., 2017;Ricciardelli et al., 2019;Parrilli et al., 2021), here we specifically investigate the molecular mechanisms through which metabolic homeostasis is maintained in a marine bacterium that (in its natural settings) is known to experience such broad seasonal temperature fluctuations. Indeed, rather than studying the consequences of a cold-shock event (a circumstance that rarely occurs in nature), we here focus on the comparison of the main cellular networks in cells growing at two different temperatures. This was performed through the integration of metabolomic and transcriptomic data unitedly with genome-scale metabolic modelling of the Antarctic bacterium Pseudoalteromonas haloplanktis TAC125 (PhTAC125). This bacterium has been isolated from an Antarctic coastal seawater sample collected in the vicinity of the T  e  r  r  e  A  d  e  l i e (66°40' S; 140° 01' E) and has received much attention in the last decade due to the interest in characterizing its cold-adaptation and nutritional adaptation strategies, as well as its biotechnological potential.
We show that different growth temperatures induce broad transcriptional changes that involve genes of many key metabolic pathways. This transcriptional rewiring, however, is scarcely reflected at the level of the core metabolism, as most key central metabolites show overlapping trends at the two tested temperatures. The obtainedomics data were used to compute the flux distributions sustaining growth at low and high temperature and this provided a mechanistic understanding of the possible adaptation strategies to temperature fluctuations.

Strains and growth conditions
PhTAC125 (Médigue et al., 2005) cells were grown in a 1.5 L GG medium (Sannino et al., 2017) (Sannino et al., 2017) in a Stirred Tank Reactor 3 L fermenter (Applikon) connected to an eZ2 Bio Controller (Applikon) at two different temperatures (0 °C -15 °C). This medium contains glutamate and gluconate as the only carbon sources. The bioreactor was equipped with the standard pH-, pO2-, level-and temperature sensors for the bioprocess monitoring. For the growths of the PhTAC125 bacterium, the preculture was centrifuged (6000 × g, 20 min, 4°C); the cells were washed twice with fresh medium, and then used to inoculate the bioreactor with a starting OD 600 of 0.2, in aerobic conditions (50% and 30% dissolved oxygen at 15 °C and 0 °C, respectively), in stirring (500 rpm at 15 °C and 250 rpm at 0 °C). The bacterial culture was carried out at 15 °C for 40 h or at 0 °C for 240 h. Each culture condition was repeated three times.
Cell growth was monitored, measuring the OD 600 about every two hours in the experiments at 15 °C, and every eight hours at 0 °C. Three different measurements were performed at each time point for each biological replicate. For intra-and extracellular metabolites analysis, the samples were taken in triplicate at five different time points during the growth in GG medium at two different temperatures (70 h, 141 h, 179 h, 190 h and 240 h at 0 °C -6.5 h, 14 h, 20 h, 25 h, and 39 h at 15 °C). For the analysis of extracellular metabolites, aliquots (1 mL) of cell cultures were harvested during the growth and centrifugated for 15 min at 13,000g at 4 °C; the supernatant was recovered, filtered (Filtropur 0.2 µm, SARSTED AG & Co. KG) and stored at -80°C. The analysis of intracellular metabolites was performed on 60 OD 600 pellets recovered during the growth by centrifuging for 20 min at 6,000 rpm at 4 °C.

Transcriptomics
For the RNA-Seq experiment, 1 OD 600 pellets were recovered during the exponential growth phase (~1/1.5 OD/mL) at 15 °C and 0 °C by centrifugation (10 min, 13,000g, 4 °C). Cell pellets were washed in RNAse-free PBS three times and stored at -80°C. Total RNA was isolated from the cells using the Direct-zol RNA Kit (Zymo Research, Irvine, CA, USA) following the manufacturer's instructions. Contaminating genomic DNA was then removed through treatment with RNAse-free DNase I (Roche, Mannheim, Germany).

RNA Sequencing
RNA concentration was measured using Quant-IT RNA Assay Kit-High Sensitivity and a Qubit Fluorometer (Life Technologies) and its quality and integrity assessed with the Agilent 4200 Tapestation System (Agilent Technologies). Indexed libraries were prepared starting from 400 ng of total RNA according to Universal Prokaryotic RNA-Seq Library Prep kit (Tecan). Final libraries were sequenced at a concentration of 1.7 pM/lane on the NextSeq 500 platform (Illumina Inc) in paired-mode 2x75bp. For each experimental condition, three biological replicates were prepared.

RNA-Seq data analysis
The twelve samples were assessed for base call quality and adapter content using fastp (Chen et al., 2018), allowing down to a mean quality threshold of 20 (i.e., probability of incorrect base call of 1 in 100) and minimum read length of 40 nucleotides. A median of 98.9% of the reads passed quality check, indicating that sequencing was carried out pristinely, and our data were biologically reliable. Salmon (Patro et al., 2017) index was built on PhTAC125's transcriptome free of tRNA and rRNA sequences, providing the entire genomic sequence as background decoy to account for possible underlying DNA contamination. Coding sequences (CDS) and genomic FASTA were retrieved on the NCBI using assembly accession GCF_000026085.1. Moreover, CDS and genomic files were concatenated with data from pMEGA plasmid (NZ_MN400773.1) and pMtBL plasmid (NZ_AJ224742.1).The transcriptome quantification step was performed using the --validateMappings flag, to ensure a sensitive selective alignment of the sequencing reads. The resulting mapping rates ranged between 56.8% and 66.21% (average of 62.47%), whereas roughly half of the reads represented ribosomal DNA sequences (discarded for downstream analyses). The integration of transcript-level abundance estimates from Salmon with the data analysis pipeline was performed using R package tximport (R Core Team, 2022, https://www.R-project.org/, R version 4.0.3). The testing of changes in the overall transcriptional output was handled through the statistical engine in R package DESeq2, with an expected proportion of false positives set at 5%. The volcano scatterplot showing statistical significance versus magnitude of change for 607 differentially expressed genes was produced using the library EnhancedVolcano from the R Bioconductor 3.14 suite.
Metabolomics 1 H NMR-based metabolomic analysis were performed on cell lysates and growth media to monitor the intracellular metabolites and the uptake and release of the extracellular metabolites, respectively, by measuring their concentration levels in samples collected at different time points during the cell growth. Medium samples were prepared in 5.00 mm NMR tubes by mixing 60 μ L of a potassium phosphate buffer (1.5 M K2HPO4, 100% (v/v) 2H2O, 10 mM sodium trimethylsilyl [2,2,3,3−2H4]propionate (TMSP), pH 7.4) and 540 μ L of each growth medium. Cell lysate samples were prepared in 5.00 NMR mm tubes by mixing 60 μ L of 2H2O and 540 μ L of samples.All the NMR spectra were recorded using a Bruker 600 MHz spectrometer (Bruker BioSpin) operating at 600.13 MHz proton Larmor frequency and equipped with a 5 mm PATXI 1H-13C-15N and 2H-decoupling probe including a z axis gradient coil, an automatic tuning-matching (ATM) and an automatic and refrigerate sample changer (SampleJet). A BTO 2000 thermocouple served for temperature stabilization at the level of approximately 0.1 K at the sample. Before measurement, samples were kept for 5 minutes inside the NMR probe head, for temperature equilibration at 300 K. 1 H NMR spectra were acquired with water peak suppression and a standard NOESY pulse sequence using 128 scans, 65536 data points, a spectral width of 12019 Hz, an acquisition time of 2.7 s, a relaxation delay of 4 s and a mixing time of 0.1 s. The raw data were multiplied by a 0.3 Hz exponential line broadening before applying Fourier transformation. Transformed spectra were automatically corrected for phase and baseline distortions. All the spectra were then calibrated to the reference signal of TMSP at δ 0.00 ppm using TopSpin 3.5 (Bruker BioSpin srl). The metabolites, whose peaks in the spectra were well resolved, were assigned and their levels analysed using a dedicated R script developed in-house. In total, 34 and 17 metabolites were identified and quantified in the cell lysate and in the growth medium spectra, respectively. The assignment was performed using an internal 1H NMR spectral library of pure organic compounds (BBIOREFCODE, Bruker BioSpin), stored reference NMR spectra of metabolites and spiking experiments. Matching between new NMR data and databases was performed using the Assure NMR software (Bruker BioSpin). The relative concentrations of the various metabolites were calculated by integrating the corresponding signals in defined spectral ranges, using in house developed R 3.0.2 scripts. Similarly, downstream data analysis was performed using R. Raw data and post-Raw data and post-processing codes are made available at https://github.com/combogenomics/MetRob015. Uptake rates of gluconate and glutamate at the two different temperatures were computed as the ratio between the average growth rate between T1 and T3 (ߤ) and the biomass yield (λ). The first was computed as: While the second was computed according to the following relationship: Biomass was obtained from OD values as described in (Giuliani et al., 2011), using 0.74 as a scaling factor for the growth at 15° and 0.66 for the growth at 0°C.

Genome-scale metabolic modelling
The genome-scale metabolic reconstruction used in this work is the one recently used in . All the simulations were run in MATLAB 2019a, using the COBRA toolbox (Heirendt et al., 2019) version 2.7.4. Metabolomic and transcriptomic data were integrated using REMI method (Pandey et al., 2019), providing TPM gene expression and relative concentration values for metabolites at timepoint T1. Both for metabolites and gene expression, we selected the top 3% as upregulated and the bottom 3% as downregulated (REMI default is 5%). Consistently with the original REMI publication, we used the two-fold change as the cut-off threshold to identify the significant gene expression and metabolite changes. All the other parameters were left as default. The predicted internal concentration of each metabolite was computed using the computeFluxSplits function implemented in the COBRA toolbox. This function computes the relative contributions of fluxes to the net production and consumption of a specific set of metabolites included in the model. Statistical tests on flux distributions were calculated using R. To estimate the difference in the activity of the main central metabolic pathways at 0° and 15°C, we computed the predicted flux of each reaction in these pathways (TCA cycle, glycolysis, Pentose Phosphate Pathway (PPP), and fatty acids metabolism) at the two different temperatures and averaged this number of reactions by the total number of reactions included in that pathway. From this latter analysis, we excluded those reactions that i) had a flux equal to zero (were inactive) in both conditions and ii) had a different sign (i.e., changed direction) at the two temperatures.

Results
We cultivated PhTAC125 cells in a bioreactor and sampled the 0 and 15°C growth curves in five different time points ( Figure 1A Figure 1A and 1B) we also performed quantified gene expression levels using RNA-Seq.

PhTAC125 metabolome is qualitatively and quantitatively robust to temperature shift
Overall, we assigned and analyzed the concentration of 34 intracellular metabolites in the two growth curves. These metabolites represented key intermediates of central metabolic pathways such as the TCA cycle, amino acids biosynthesis, glycolysis, Pentose Phosphate Pathway (PPP), and nucleic acids biosynthesis. Our NMR metabolomic approach was not able to discriminate between the oxidized and reduced forms of NAD and NADPX so we generally refer to NADX and NADPX in the rest of the manuscript. First, we monitored the overall trends of their intracellular concentrations with respect to the beginning of the growth experiment ( Figure 2A) by computing the Pearson product moment for each metabolite across the five time points and the relative statistical support (Spearman correlation, p-value<0.05, Figure 2B). Strikingly, for about 90% of the analyzed metabolites we found a positive correlation between the variation of their concentrations during the two separate growth experiments (at 0 and 15 °C). These positive values ranged from 0.21 in the case of aspartate to 0.96 in the case of malate. Ethanol and thymidine displayed a Pearson Product-Moment correlation (PPM) quite close to zero (-0.16 and -0.1, respectively), thus showing no correlation in the two experiments. Finally, the concentrations of PEP and ribose were negatively correlated in the two growth experiments, with PPM of -0.9 and -0.57, respectively. The negative correlation of PEP was supported statistically (Spearman correlation p-value = 0.03692).
To qualitatively assess the change in concentration of each metabolite against each other, we computed the all-against-all correlation among the metabolites identified in our study. The results of this analysis are reported in Figure 2D. We observed two main blocks of metabolites, one embedding amino acids (with the exception of glutamate and glutamine) and one mostly including intermediates of the main central metabolic pathways (TCA, PPP, nucleic acids biosynthesis). Each of these clusters was characterized by a strong (and statistically supported) correlation among the representatives of the same cluster and an equally relevant anti-correlation with the members of the other cluster. Overall, the intracellular concentration of amino acids was shown to increase over time, whereas the level of, for example, TCA and PPP intermediates was shown to decrease throughout the growth experiment. To assess whether the two growth temperatures had a role in determining the observed trends in metabolite concentrations, we computed an all-against all-correlation for each metabolite at 0° and 15°C separately and then evaluated whether these two correlation matrices differed significantly. The t-tests performed supported no significant differences between the means of these two comparisons (p-value = 0.74 and 0.89, respectively) suggesting that, despite the large difference in growth temperatures, the overall dynamics of intracellular core metabolites was maintained, pointing to an apparent structural metabolic robustness to growth temperature. This robustness was also conserved quantitatively, as the average concentration of each metabolite was maintained similar across the two different growth experiments ( Figure 2C). Indeed, with the exception of 6 metabolites out of 34 (Acetate, Glutamine, Glucose, PEP, NADX and NADPX, t-test, p-value < 0,05), no other average metabolite concentration showed a significant difference at 0 and 15°C ( Figure 2C and Figure S1). A similar scenario was observed for the pool of (17) quantified extracellular metabolites ( Figure 3A). Consistently with growth conditions, glutamate and gluconate concentrations decreased over time until their exhaustion from the growth medium ( Figure 3A), a trend i) that indicated the almost simultaneous consumption of these two nutrients and ii) that was observed at both 0 and 15°C. Out of the 17 extracellular metabolites that were quantified in our experiments, 16 showed a positive correlation between their trends at 0 and 15°C ( Figure 3B). In particular, 10 of them displayed a PPM above 0.5. The only metabolite whose trend differed in the two experiments was 2-oxoglutarate (PPM = -0.4). This situation was mirrored at the quantitative level ( Figure 3C). Indeed, none of the extracellular metabolites showed a statistically significant difference between 0° and 15°C growth experiments. Finally, metabolomic data allowed us to evaluate the uptake kinetics of glutamate and gluconate in the tested conditions ( Figure 1C and D) revealing that the two carbon sources were taken up almost simultaneously in the two experiments and were depleted/exhausted around T4.

The 0 and 15°C transcriptomes of PhTAC125
We next characterized the transcriptomes of PhTAC125 at 0 and 15°C growth temperatures. We sampled both curves during exponential growth (sample points "1" in Figure 1A and B) and sequenced the transcriptome using RNA-Seq technology (see Material and Methods section). The main features of transcriptomic data are shown in Table S1 and Figures S2 and S3. Differentially expressed genes (DEGs) were identified using a log 2 fold change of 1 (or -1) and an adjusted p-value of 0.05 as thresholds. Overall, we identified 607 differentially expressed genes in the comparison between growth at 15° and 0° C, with an almost identical amount of up-and down-regulated genes at 15°C (304 and 303 genes, respectively). The outliers of up-and down-regulated genes are shown in Figure 4A and described in Table 1. Remarkably, out of 607 DEGs, 359 were metabolic protein-coding genes (roughly 59%).
Here we first describe specific features of outlier DEGs and then focus our attention on broader functional categories that embedded the highest fraction of differentially expressed genes. Among the top-five downregulated genes, two lacked a clear functional annotation (PSHA_RS12380 and PSHA_RS08940) and further experiments are required to understand their role in the adaptation to growth at warmer temperatures. Two of them (PSHA_RS07950 and PSHA_RS06965) were annotated as TonB-dependent receptors, suggesting their involvement in the uptake and transport of large substrates, possibly siderophores complexes and/or vitamins. The remaining gene (PSHA_RS11830) is annotated as NAD-dependent succinate-semialdehyde dehydrogenase, responsible for the conversion of succinate-semialdehyde to succinate. The downregulation of this gene is in line with the lower intracellular concentration of succinate at 15°C, as discussed later in the text. Strikingly, among the over-expressed genes during growth at 15° C, we found a prevalence of cold-shock related proteins (4 out of five) and a ribonuclease R encoding gene (PSHA_RS14645).
Although no functional studies have been done on cold shock responses and coldshock proteins (CSPs) from psychrophilic bacteria, similarities with the CSPs that are produced in mesophiles have been observed. In particular, the over-expressed CSPs in PhTAC125 share an amino acid sequence identity ranging between 63.8% and 65.6% with the mesophilic cspE of E. coli. Furthermore, these proteins contain highly conserved RNA-binding motifs, RNP1 (K-G-F-G-F-I) and RNP2 (V-F-V-H-F) (Landsman, 1992;Yamanaka et al., 1998), indicated by a black box in Figure S4.
PSHA_RS14630, PSHA_RS14635 PSHA_RS14640 and PSHA_RS16715 contain one highly conserved nucleic acid-binding domain called cold-shock domain (CSD, http://pfam.xfam.org/family/PF00313), which is annotated as RNA chaperone/antiterminator. Moreover, the NCBI identifies both genes PSHA_RS14635 and PSHA_RS14640 as encoding for the same protein (WP_011329604.1), while the others encode for proteins with a different identifier but virtually same sequence. Cold shock proteins would normally counteract the deleterious effects of temperature drop, enabling the cells to grow at low temperatures (Keto-Timonen et al., 2016). Jiang et al proposed the role of cspA as an RNA chaperone capable of melting the RNA secondary structure (Jiang et al., 1997), thereby enhancing translation of mRNAs at low temperatures. However, counterintuitively with their given name, not all members of the CSP family are cold-inducible and their expression is activated upon different stresses (Caballero et al., 2018). For this reason, CSPs might be required for bacterial adaptation to environmental changes. The fifth upregulated protein, Ribonuclease R, belongs to the RNR family. In Escherichia coli, the Ribonuclease R consists of a central nuclease domain, two cold-shock (CSD) domains near the N-terminal region of the protein, an S1 domain and a highly basic C-terminal region (Phadtare and Severinov, 2010). Cairrão et al showed that the rnr gene is co-transcribed with flanking genes as an operon induced under cold shock in E. coli (locus tag b4179) (Cairrão et al., 2003). There is an important analogy here since, although PhTAC125's operon map isn't available, its top 4-upregulated genes are mapped consecutively on the genome, with coding sense on the same strand (+), spaced by an average distance of 216 nucleotides. The homology to a variety of domains involved in stress response, together with our observation of a strong upregulation during growth at 15°C may indicate that PhTAC125 strives to mediate the elimination of detrimental secondary structures and a temperature rise promotes the expression of enzymes that are required for the correct processing of rRNA precursors. Extending the analysis to a broader functional level, we identified significant differences in 12 COG categories. Specifically, intracellular trafficking, secretion and vesicular transport, signal transduction mechanisms, cell motility, cell wall/membrane/envelope biogenesis, replication, recombination and repair and cell cycle control, cell division, chromosome partitioning functional categories were found to be over-represented among up-regulated genes during growth at 15°C. This matches the results from recent works on the characterization of the evolutionary pathways responsible for thermal adaptation and the overall notion that one of the main effects of temperature increases on cell physiology consists in the disruption of membrane integrity caused by increased fluidity (Angilletta Jr., 2009;Toll-Riera et al., 2022). Overexpression of cell wall and membrane biogenesis-related genes might help to overcome this feature associated with growth at higher temperatures. Conversely, genes involved in post-translational modification, protein turnover, chaperones, translation, ribosomal structure and biogenesis, lipids, nucleic acids and amino acids transport and metabolism as well as energy production functional categories were over-represented among down-regulated genes during growth at 15°C ( Figure 4B). We believe that this analysis depicts the system-level adaption of bacterial life to diverse growth temperatures and is in line both with the physiological features observed in this work and with previously obtained data on cold/warm adaptation. Indeed, the significant over-expression of genes broadly related to cell replication (COG categories L, M and D) is in line with the higher growth rate observed at 15°C in respect to 0°C, over the entire growth period ( Figure 1A and B) and in the specific time interval where RNA was sampled. Further, the over-expression of cell motility at warmer temperatures is in perfect agreement with previous assays on PhTAC125 motility that had shown a reduced swimming capability of this strain at 0 °C (Ricciardelli et al., 2019). Surprisingly, the up-regulation of cell-growth related processes is counteracted by a general down-regulation of cellular metabolism. Key processes of PhTAC125 core metabolism such as amino acid, lipid, and nucleotide metabolism together with energy production/conversion were found to be down-regulated ( Figure 4C). As mentioned previously, among 607 DEGs, 359 (59.14%) were metabolic genes. Zooming in at the level of the single pathways ( Figure 4C) revealed that the downregulation affected most of the key central pathways of PhTAC125 metabolic network (including, for example, glycolysis, TCA cycle, PPP and amino acids biosynthesis genes). In these pathways the number of downregulated genes strongly outpaced that of up-regulated genes ( Figure 4C). Thus, despite consistent intracellular and extracellular metabolic profiles between the two tested growth temperatures, the underlying expression of metabolic genes showed an opposite trend, with a remarkable number of DEGs in the two growth conditions. For example, the intracellular pools of key TCA intermediates (fumarate and malate, Figure 2C) were shown to be similar between 0 and 15 °C, despite more than 50% of the TCA genes were shown to be differentially expressed in the two conditions ( Figure 4C). Similarly, the intracellular concentration of amino acids was shown to be comparable between the two conditions ( Figure 2C) whereas the level of expression of their corresponding biosynthetic genes was shown to be significantly different ( Figure  4C). Accordingly, we hypothesize that the transcriptional network of PhTAC125 provides a buffering mechanism through which metabolic homeostasis is maintained, at least at the level of the central metabolism. In the next session we will combine metabolomic and transcriptomic data with a genome-scale metabolic reconstruction of PhTAC125 to unravel the intimate mechanisms through which this balance is achieved.

Genome-scale modelling of growth temperature adaptation
A genome-scale metabolic reconstruction exists for the strain PhTAC125 (Fondi et al., 2015). We thus exploited this resource to get a mechanistic interpretation of metabolic homeostasis. First, we checked whether the model was able to represent the experimentally determined phenotypes. We thus constrained the PhTAC125 genomescale reconstruction with uptake rates of glutamate and gluconate measured at 0 and 15° degrees and run an FBA simulation to predict the cellular growth rates. Average uptake rates of glutamate and gluconate were respectively computed between T1 and T3 for both growth curves (that is after 14 and 141 hours for 15 and 0°C, see Material and Methods), and resulted in 0.10 and 0.08 mmol/gCDW* -1 at 0°C and 0.47 and 0.62 mmol/gCDW*h -1 at 15°C. As shown in Figure 5A, at 0°C the results of this simulation (0.021h-1) were in line with the measured growth rates (0.027h -1 ). Conversely, at 15°C a discrepancy was observed between the experimental growth rate and the one predicted in silico, 0.27 vs 0.13 h -1 , respectively. Although there could be many possible explanations for this, we hypothesize that constraining the model only based on the uptake rates of the provided nutrients may not be enough to correctly represent the metabolic phenotype of the cells exposed to different temperatures. For this reason, and to provide a more realistic picture of how metabolic homeostasis is maintained at the level of the central metabolism, we combined metabolomic and transcriptomic data with the genome-scale modelling of PhTAC125 metabolism. We decided to focus on T1 (beginning of the exponential growth phase) as this time point should better resemble the cellular physiological state in which FBA assumptions hold the most (i.e. metabolic steady state). Then, we computed the log 2 fold change (log 2 FC) of central intracellular metabolites in the two conditions ( Figure 5C). Most (76%) of the metabolites showed a | log 2 FC | lower than 1, confirming an overall robustness of the central metabolism to growth temperature. Only 8 metabolites displayed a | log 2 FC | lower/greater than 1 in the contrast between 15° and 0° C metabolomic data and for only three of them we obtained a statistical support (p-value<0.05, Figure 5C). More specifically, NADX and Gluconate were the two metabolites showing a log 2 FC > 1 (thus being more abundant at 15°) whereas PEP, Succinate, UDP-GlcNAc, Thymidine, acetate and X1-methylnicotinamide displayed a log2FC < -1 (thus being more abundant at 0°).
We then constrained the model using all available experimental data obtained in this work using REMI (Pandey et al., 2019) (see Material and Methods). Briefly, this approach allows translating gene expression and metabolite abundance data resulting from a "perturbation" experiment into differential flux distributions among the two resulting conditions. In our case, we used data (metabolomes and transcriptomes sampled at the same time point) from 15° and 0° C growths to analyze the systemic metabolic alteration(s) occurring in this pair of conditions. FC and TPM values ( Figure  5B and C) were used to constrain the model with metabolic and expression data, respectively. Further, the model was also constrained by setting the boundaries of uptake reactions to represent the actual medium used for the experiment (Schatz medium with glutamate and gluconate as the sole carbon sources, see Materials and Methods). The outcome of data integration into the metabolic reconstruction is summarized by REMI through the computation of the theoretical maximum consistency score (TMCS) and the maximum consistency score (MCS). The first indicates the number of available omics data (for metabolites and reactions) whereas the latter represents the number of those constraints that are consistent with fluxes and could be integrated into REMI models. As a result, the MCS is always equal to or smaller than the TMCS. In other words, MCS is the largest fraction of available data (metabolomics and transcriptomics) that could be incorporated into an FBA model from a given set of constraints (the abundance of metabolites and transcripts) while ensuring that the model still achieves the required metabolic functionalities and remains feasible. TCMS indicates the number of genes and metabolites with available relative abundance values that can potentially (either because above the specified threshold or their actual inclusion in the metabolic reconstruction) be integrated into the model. The proportion between TMCS and MCS obtained in this work is comparable to that from other studies where REMI was used for the same purposes (Pandey et al., 2019) ( Figure 5D). Out of data integration and FBA simulations, we obtained two distinct flux distributions, i.e. the flux distribution resembling growth at 0°C and the one theoretically accounting for the growth at 15°C. The differences between the two will represent the most likely metabolic alterations in response to growth temperature and, consequently, will highlight those pathways/reactions that contribute the most to maintaining the observed metabolic homeostasis. First, we checked whether our simulations were accounting for the actual differences between 0° and 15 °C intracellular metabolomes. We thus computed a matching coefficient (Simple Matching Coefficient, SMC) between the predicted and measured log 2 FC of each of the 34 internal metabolites, accounting for how many times the model correctly predicted the increase (or decrease) of its internal concentration. Overall, we found an SMC of 71% between simulated and measured metabolic data, revealing that the model is capable of accounting for most of the central metabolome rewiring in respect to growth temperature. Then we focused on those metabolites that showed a marked change between growth at 15 and 0°C (i.e. for which | log 2 FC | > 1, Figure 5C). Figure 5F shows that, except for three metabolites (UDP-GlcNAc, Thymidine and X1-methylnicotinamide) for which the model does not predict any difference between the two growth conditions, for the remaining 5 metabolites the model correctly predicts a higher (NADX and gluconate) or lower (acetate, PEP and succinate) internal production at 15°. Overall, we found a significant positive correlation of 0.72 between the measured and the predicted FC of internal concentrations of the 8 metabolites for which | log 2 FC | > 1 ( Figure 5F) (Spearman correlation, p-value= 0.04).
Next, we examined at what extent the central, interconnected pathways showed significantly altered flux distributions. While no clear signal could be identified for glycolysis and TCA cycle (mirroring what we observed with DEGs analysis), both PPP and fatty acids metabolism showed significant differences in their fluxes between 0° and 15° growth simulations. For what concerns fatty acids metabolism, fluxes representative of growth at 0°C were significantly higher (Kolmogorov-Smirnov test, pvalue = 8.96 -10 ) than those resembling metabolism at 15°C ( Figure 5F). The opposite was observed for PPP simulated fluxes, which displayed, on average, significantly lower values at 0° vs 15°C (Kolmogorov-Smirnov test, p-value = 0.0030). A sustained activity of PPP at 15° is in line with the increased internal concentration of gluconate and ribose (two key PPP intermediates) at this temperature ( Figure 5C). Similarly, an increase in fatty acids biosynthesis (and fatty acids metabolism in general) is in line with previous simulations and experiments concerning the involvement of this process in the adaptation to growth at (relatively) low temperatures (Fondi et al., 2015). We then focused on the analysis of simulated metabolic fluxes around the metabolite that showed the higher degree of variation, i.e. PEP. We asked which metabolic rewiring could lead to an increased production of this metabolite at 0°C. The analysis of fluxes revealed an increased activity at 0°C with respect to 15°C of both the reaction converting L-Malate to 2-Oxaloacetate ((S)-malate:NAD+ oxidoreductase) and the reaction converting 2-Oxaloacetate to PEP (phosphate:oxaloacetate carboxy-lyase). Also, ATP:pyruvate,H 2 O phosphotransferase reaction displayed an increase in flux at 0°. Overall, this would allow the partial redirection of key TCA cycle intermediates to the production of PEP. At the same time, we also recorded an increased flux at 0°C in the reaction redirecting Acetyl-CoA to fatty acids-biosynthesis. Conversely, reactions leading to the production of citrate from OAA and Acetyl-CoA resulted to be less active at 0°C. Overall, our experimentally constrained simulations seem to suggest a working model for the metabolic adaptation to growth at different temperatures. Considering growth at 0°C, an increased level of intracellular PEP observed experimentally might be functional to its conversion to Acetyl-CoA and its consequent tunneling into fatty acids metabolism ( Figure 5G) and might be the outcome of fluxes redirection from the TCA cycle to PEP production, rather than from glycolysis or PPP. Consistent with this idea is the overall increase of key TCA intermediates ( Figure 5A) measured during growth at 0°C with respect to growth at 15°C that suggests an overall increased activity of this central pathway in the 0°C growth condition.

Discussion
In this work, we have studied how a cold-adapted bacterium rewires its central metabolism when growing at two distinct temperatures that overall resemble a seasonal shift. This was done by characterizing a pool of 34 intracellular and 17 extracellular central metabolites during 5 different time points of its growth curves (i.e. at 0°C and 15°C) and by evaluating gene expression during the initial stages of its exponential phases. To our surprise, the top five upregulated genes at 15°C contain a cold shock domain, which is usually recruited to counteract the deleterious effects of temperature drop. A possible explanation for the apparent paradox where CspA is activated not only following cold stress but also under non-stress and other stress conditions which entail a down-regulation of bulk gene expression and protein synthesis is presented in (Brandi and Pon, 2012). CSPs have been found in almost all types of bacteria and are mainly induced after a rapid temperature downshift to regulate the adaptation to cold stress but are also present under normal conditions to regulate other biological functions (Horn et al., 2007). For example, in E. coli, only four (cspA, cspB, cspG, and cspI) of the nine CSP genes are cold-induced (Xia et al., 2001). Two of them, cspE and cspC, are constitutively expressed at physiological temperatures, and act as 'housekeeping RNA chaperones' to modulate the global gene expression (Phadtare and Inouye, 2001;Phadtare et al., 2006). Furthermore, these CSPs are also involved in the transcription anti-termination mechanism, which is based upon preventing the formation of secondary structures on the nascent mRNA (Bae et al., 2000;Phadtare et al., 2002;Phadtare and Severinov, 2005). In PhTAC125, four cold shock-like proteins cspE are upregulated during growth at 15 °C. This indicates that these CSP proteins may not play a role in cold adaptation in PhTAC125, but analogously to E. coli, they may act as chaperones by destabilizing secondary structures in target RNA at high temperature so that the single-stranded state of target RNA is maintained. This may then enable efficient transcription and translation.
From a broader perspective, the comparison between the 0° and 15°C transcriptomes produced more than 600 differentially expressed genes, assigned to a dozen of different functional categories. On the other hand, the intracellular concentration of nearly 90% of the analyzed metabolites correlates positively in the two conditions. We interpret this as the capability of PhTAC125 regulatory network to buffer the temperature shift to produce strikingly similar metabolic phenotypes, despite the temperature gap. This is also reflected in the conservation of the overall metabolic network structure (i.e., the presence of correlated clusters of metabolites) ( Figures 2C  and 3C). Remarkably, this conservation involved metabolites from many different pathways such as amino acids metabolism, TCA cycle, PPP and glycolysis, suggesting that such robustness is propagated at the level of the entire metabolism and is not restricted to specific pathways. Studies on the metabolic response to growth at different temperatures in other microorganisms have generally shown larger variability among intracellular pools of metabolites (Guan et al., 2018) and/or the tendency to activate genetic mechanisms that shut down metabolism under longer-term high-temperature stress (Kim et al., 2020). This is not the case for PhTAC125 that, through a global transcriptional buffering maintains very consistent trends of its intracellular and extracellular metabolomes. This mirrors, for example, the response of E. coli to genetic and environmental perturbations and its capability to maintain metabolite levels stable, reflecting the rerouting of fluxes in the metabolic network (Ishii et al., 2007). A few metabolites stood out in this conservation of metabolite pools, including PEP and ribose (showing opposite trends at 0 and 15°C) and ethanol and thymidine (showing no correlation between the two experiments). Also, NADX and NADPX levels differed in the two conditions, as shown in Figure 2A.
To unravel this complex interplay and produce a mechanistic interpretation of this adaptation, -omics data (together with data on carbon sources uptake rates) were used to constrain a genome-scale metabolic reconstruction and derive the most-likely metabolic phenotypes at these two temperatures. This approach provided evidence that most of the metabolic adaptation is probably played at the level of the phosphoenolpyruvate (PEP)-pyruvate-oxaloacetate node and specifically involves the increase of TCA fluxes and their redirection to PEP production. Consistent with this idea is the overall reduction of key TCA intermediates at 15°C ( Figure 5A) and the observation that, upon growth on non-glycolytic compounds (e.g. glutamate and gluconate) the cycle intermediates malate or oxaloacetate must be converted to pyruvate and PEP for the synthesis of glycolytic intermediates. In our case, however, the availability of an increased pool of PEP would rather allow its conversion to acetyl-CoA and its consequent tunneling into fatty acids metabolism ( Figure 5G). Also the overexpression of genes involved in valine, leucine and isoleucine degradation ( Figure  4C) at 0°C can be associated with the higher activity of fatty acids metabolism, as the degradation of these two amino acids leads to the production of acetyl-CoA (Massey et al., 1976). Consistently with this observation, our model predicts an increased activity of fatty acids metabolism at 0° ( Figure 5G). The importance of fatty acids metabolism in growth at low temperatures is largely known. For example, the fluidity of the rigidified membrane (Los and Murata, 2004) imposed by low temperatures can be restored through the modulation of the i) saturated and unsaturated fatty acids, ii) fatty acid chain length, iii) the proportion of cis to trans fatty acids thus, ultimately, through the modulation of fatty acids metabolism (as observed in this work). This also validates and extends previous findings on general cold adaptation strategies (Graumann and Marahiel, 1999;Ting et al., 2010;Fondi et al., 2015), pinpointing the role of specific metabolic reactions in this process.
We also found an increased PPP activity at 15° with respect to 0°C ( Figure 5G). Increased PPP activity signifies increased NADPH production and is linked with many biosynthetic processes. Interestingly, high activity of the PPP pathway has been described as an adaptative mechanism (in soil bacteria) to temperature stress (Bore et al., 2017), in particular by linking the requirement of sugar units for biofilm formation to the overall activity in this metabolic pathway. We know for a fact from the work of (Ricciardelli et al., 2019) that PhTAC125 displays strong propensity to form more biofilm at 15°C vs. 0°C when grown in the same minimal medium used in this work. Also, we found a remarkable fraction of overexpressed biofilm-related genes ( Figure  4C), reinforcing the idea of a temperature-mediated switch to this kind of lifestyle. Thus, we infer that a similar mechanism is at play in this marine bacterium where the redirection of fluxes in PPP at 15°C would be the metabolic basis for the formation of cellular aggregates at higher temperatures. Moreover, the increased PPP activity can respond to other cell requirements linked to the DNA replication recombination and repair, since the ribose produced by PPP is necessary for these cellular activities that, as suggested by transcriptomic analysis (Figure 5B), result more active at 15°C (Table  S2). Very recently, Macarena Toll-Riera et al. have published (Toll-Riera et al., 2022) a work in which PhTAC125 was evolved at increasingly higher temperatures to study the evolutionary potential of upper thermal tolerance and characterize the genomic basis of temperature adaptation. Interestingly, the authors found almost no metabolic genes among those that apparently provided an increased fitness at higher temperatures. Indeed, clones that were selected for their improved growth at higher temperatures (up to 30°C) had mainly mutations in proteases-encoding genes or in genes involved in chromosome copy number reduction, energy production and conversion and cell wall biosynthesis. Thus, apparently, the metabolism of PhTAC125 was not the primary target of mutations that led to an increased fitness in the tested conditions. This, in turn, may indicate the presence of an already optimized and plastic metabolism that allows growth in a broad range of temperatures, in line with the metabolic robustness that we have characterized in this work.

Conclusions
The aim of our experiments was to evaluate the effect of growth temperatures to cellular homeostasis, particularly at the metabolic level. We have shown that a coldadapted bacterium expresses two very similar metabolic phenotypes in response to two widely different temperatures, by adjusting the expression of key enzymes and fine-tuning the intracellular concentration of key intermediates. This indicates that i) marine microbes can successfully adapt to broad (seasonal) temperature changes by evolving regulatory circuits that buffer such variations and maintain the underlying metabolic network robust to fluctuations and that ii) integrated, multi-omics approaches are essential to comprehend adaptations to environmental fluctuations at the whole cell level.