Hepatic Transcriptomic Responsiveness of Polar Cod, Boreogadus saida, to Ocean Acidification and Warming

Abstract

Background: This study was part of a larger comprehensive project (BIOACID) addressing the physiological resilience of Polar cod, Boreogadus saida, to ocean acidification and global warming and aimed to unravel underlying molecular mechanisms of the observed physiological responses. Methods: Fish were acclimated long-term to three CO₂ concentrations comprising control conditions (390 ppm) and two projected climate scenarios (780 ppm and 1170 ppm). Each CO₂ treatment was combined with four temperatures: 0, 3, 6, and 8 °C. Here, we focused on the hepatic transcriptomic profiles from these previously physiologically characterized fish. Results: Generally, we did not detect signs of a classical stress response. Consistent with functional observations, warming induced much stronger molecular responses compared to elevated PCO₂, but an interaction between both factors existed to some extent. Gene ontology analysis revealed a strong response in lipid, amino acid, and protein metabolism. With increasing temperature, we observed a shift away from lipid metabolism, while carbohydrate metabolic pathways remained stable. Conclusions: Although we found Polar cod to be quite resilient to ocean acidification, temperature will remain a critical parameter for this valuable Arctic keystone species, and the question remains as to whether the observed acclimation strategies can be implemented in its natural habitat, especially when food supply is limited.

1. Introduction

The Earth’s climate has changed dramatically in recent decades. Numerous studies have already identified serious changes, e.g., a substantial increase in greenhouse gas emissions such as carbon dioxide, which causes atmospheric and sea surface warming, resulting in melting of glaciers and loss of (sea)ice (e.g., [1,2,3,4]); and ocean acidification (OA, the process of rising CO₂ level and consequent dropping pH, e.g., [5,6,7,8]). The Arctic is supposed to be the region most affected by the ongoing climate change. A recent analysis revealed a warming rate four times faster than the global average in the last four decades [9]. Combined with the function of the world’s oceans as a carbon sink (which took up 2.8  ±  0.4 GtC��yr⁻¹ during the last decade [10]), the marine environment in the Arctic appears to be an extremely endangered ecosystem. It is, therefore, of utmost importance to understand the responses to environmental changes of species from this region, identify thresholds and critical tipping points, as well as draw conclusions how the organisms might already have adapted to the current changes.

Based on their ability to regulate their acid–base household, it was assumed that fish are very tolerant towards OA, but further studies have shown that OA strongly affects young life stages and spawners in particular [11,12,13], as well as the population dynamics of fish stocks [14,15]. Additionally, the high partial pressure of CO₂ (PCO₂) may affect the tolerance to other stressors (reviewed by [16,17]), and is even conjectured to constrict the thermal tolerance window in marine ectotherms [18,19,20]. However, the response to temperature is much more pronounced in fish, and can be detected on physiological, metabolic, cellular, and molecular levels. Species distribution, especially, is bound to a specific temperature range [21,22,23,24]. In contrast to species that occur over a wider latitudinal range, the adaptive potential to temperature of high-latitude species appears to be comparably low, resulting in narrower windows of thermal tolerance (see [25,26,27]). These species generally experience lower than optimal temperatures in contrast to species that inhabit lower latitudes, which are more likely to experience average temperatures closer to their upper thermal limits [22,28]. Thus, progressive warming has already led to a poleward species shift [29,30,31,32], resulting in competition between invading and endemic species and, in many cases, increasing the predation pressure on the latter. In particular, the northward migration of Atlantic cod, Gadus morhua, haddock, Melanogrammus aeglefinus, wolffishes, Anarhichas spp., redfishes, Sebastes spp., and Atlantic mackerel, Scomber scombrus, is threatening Polar cod stocks [33,34,35,36,37,38].

Polar cod is a relatively small fish of 40 cm in maximum length, belonging to the family of Gadidae [39,40]. The core thermal habitat of Polar cod is life-stage dependent, and adults usually inhabit regions between −1.5 and 2 °C [41,42]. It has an Arctic circumpolar distribution in waters with and without drifting sea ice [43,44], feeds on small invertebrates [45,46], and is itself an important prey species for various taxa [47,48]. Polar cod thus forms the link between lower trophic levels and top predators and is, thus, considered a key species in the Arctic ecosystem (see also [34],reviewed by [46,49,50]). Due to its role as a keystone species, the interest in Polar cod’s ecology, physiology, and resilience to climate change has been the focus of several studies (e.g., [34,46,51,52,53,54,55]), including the multi-stressor BIOACID project on the combined effects of ocean acidification and warming on marine biota (www.bioacid.de, accessed on 15 May 2024), in which our study was involved [56,57,58,59,60,61]. In this comprehensive project, a variety of physiological parameters such as the energy budget, growth, and exercise [56,57]; metabolic background of performance in terms of cardiac mitochondrial function [58]; as well as neurochemical and behavioral shifts [60] of Polar cod were investigated in response to ocean warming and acidification. Unlike the physiological studies, the present transcriptomic study includes animals treated at intermediate PCO₂ (780 ppm) to determine possible molecular fine-tuning due to progressing OA. For the majority of the tested physiological parameters, the overall temperature effect was more pronounced, and progressing PCO₂ led predominantly to only slight trends or remained under the detection limit. The experiments indicated that the optimum temperature for Polar cod is around 6 °C, at which both growth performance and aerobic scope are maximal [56,57], which is further corroborated by the optimum temperature for ATP production in the range of 3–6 °C [58]. This indicated a sufficient oxygen supply capacity even under activity at 6 °C (cf. [57,62] under review). Thermal limitations were detected at the highest applied experimental temperature, indicating a long-term upper thermal limit at 8–10 °C. At this temperature, oxygen demand (here measured as standard metabolic rate) and mortality increased and growth decreased, accompanied by a reduced feed conversion efficiency [56]. Cardiac mitochondrial efficiency decreased as well, indicated by an increased proton leak, leading to decreased ATP production efficiency [58]. This was in line with significant temperature effects on the neurochemical composition, which particularly evoked a decrease in amino acid and osmolyte metabolism, underlining the energetic limitations of Polar cod at high temperatures [61] and potential shifts in energy allocation due to a disrupted balance between energy production and demand. Behavioral studies revealed a hypercapnia-induced, temperature-independent reduction in laterality, suggesting reduced fitness of Polar cod under future CO₂ conditions [60]. This also became visible in their impaired swimming performance under elevated PCO₂ [57].

The basis of all the above-mentioned studies of the BIOACID project was a large cross-factorial acclimation experiment in which fish were acclimated long-term to three different CO₂ concentrations: 390 ppm (control), a moderate increase (780 ppm), and a high PCO₂ corresponding to the business-as-usual greenhouse gas emission scenario (RCP8.5, 1170 ppm, [63]). Each PCO₂ level was combined with four temperatures—0, 3, 6, and 8 °C—to simulate different warming scenarios and to test for physiological acclimation thresholds. In this study, we focused on the hepatic transcriptomic profiles from these fish to identify relevant molecular mechanisms underlying the previously reported physiological responses.

2. Materials and Methods

2.1. Animal Maintenance

The fish were caught in January 2013 by RV Helmer Hanssen and kept at the Havbruksstasjonen i Tromsø AS (HiT) at 3.3–3.8 °C for three months before being transferred to the aquaria of the Alfred-Wegener-Institute (AWI) in Bremerhaven. Fish were then maintained in flow-through tanks at normoxia and an ambient temperature of 0.0 ± 0.5 °C. Four weeks prior the start of the experiments, the temperature was increased to 5 °C. During this time of preconditioning, the fish were fed daily with high-protein food pellets (Amber Neptun, 5 mm, Skretting AS, Stavanger, Norway, [56]).

2.2. Experimental Setup: BIOACID Incubations

Detailed information on the experimental conditions and incubation setup can be found in Kunz et al. (2016) [56] and Leo et al. (2017) [58]. The full raw water chemistry data of the incubation can be found in Schmidt et al. (2016) [64].

Briefly, four different temperatures were adjusted with a maximum rate of 2 °C/d until the incubation temperatures of 0, 3, 6, and 8 °C were reached. Each temperature was combined with one of three CO₂ levels: 390 (control scenario, current habitat conditions (2013)), a moderate increase of 780 [RCP6.0 according to, [63]] and a high increase of 1170 ppm CO₂ (business-as-usual greenhouse gas emission scenario, RCP8.5, according to [63]). Due to a limited number of available specimens, the treatment of 0 °C and 780 ppm was omitted. For each combined treatment, 12 single aquaria (approx. 24 L each) were used, and a dark–light cycle of 12: 12 h was maintained. The light was dimmed during the day, except when feeding and removing food and while the water condition was measured. The respective CO₂ levels were pre-adjusted in a header tank containing ~200 L of seawater, supplying the individual aquaria. A mass flow controller (4- and 6-channel MFC system, HTK, Hamburg, Germany) was used to set the desired PCO₂ by mixing almost CO₂-free pressurized air with pure CO₂. After four months of incubation, a set of physiological experiments was carried out [56,57,58,60,61] for a subset of animals, while the remaining animals were sampled directly for molecular and biochemical analyses (sampling round 1: procedure as described below). After 2 weeks of recovery, the experimental fish were anesthetized with 0.2 g·L⁻¹ tricaine methane sulphonate (MS-222) (sampling round 2). After the fish were unconscious, they were weighed, standard and total length were determined, a blood sample was taken, and hematocrit was determined. Afterwards, fish were killed by a cut though the spine; all organs were sampled; and the liver, gonad, stomach, and gutted weights were determined. All tissue samples were immediately snap-frozen in liquid nitrogen and stored at −80 °C. During optimization of the microarray (see Section 2.5), no differences between the sampling rounds could be detected. To exclude a bias based on sex differences, only males were used from both sampling rounds for the microarray experiment (n = 5 per combined treatment).

2.3. Whole Animal Energetics

2.3.1. Growth Experiment

During the incubation period, the growth of the fish was monitored by measuring individual fish length (total length, to the nearest mm below) and weight (wet weight, to the nearest 0.1 g) at three time points. The measurements took place on the third day after feeding. Detailed information about the procedure can be found in Kunz et al. (2016) [56].

Briefly, the first measurement took place prior to the transfer into the experimental setup (day 0); the second after half of the incubation time had passed (day 65); and the third and last at the end of the incubation (day 130), right before the tissue sampling. Before the measurements, the fish were anesthetized with a light dose of MS-222 (0.06 g L⁻¹).

The fish were fed ad libitum once every four days with a with a pre-weighed amount of high-protein feed pellets (relative nutrient composition: 54.7% protein, 19.1% fat, 8.0% carbohydrates, and 8.5% moisture; Amber Neptun, 5 mm, Skretting AS, Norway), which were soaked in a specific amount of sea water 24 h prior feeding for better feed acceptance. The aquaria were cleaned daily before feeding.

Specific Growth Rate (SGR)

The specific growth rate (SGR) was calculated in percent according to Jobling (1988) [65] and defined as the rate of body weight increase per day:

$$\mathrm{S}\mathrm{G}\mathrm{R}=\frac{(\ln \left({\mathrm{W}}_2\right)-\ln \left({\mathrm{W}}_1\right))}{{\mathrm{t}}_2{-\mathrm{t}}_1} 100$$

with W₂ and W₁ being the individual weights in grams at day one and day two measurement, respectively, and t₂ − t₁ being the time between weighings in days.

2.3.2. Physiological Condition before Sampling

In this study, Fulton’s condition index (FCI, [66,67]), the gonadosomatic index (GSI, [68]), and the hepatosomatic index (HIS, [69,70]) were calculated in order to characterize each individual’s fitness and to correlate the effects of the treatments on the condition and growth.

The FCI is often used to estimate the whole-body condition and indicates changes in the food reserves stored in muscle [71]. According to Fulton (1904, 1911) [66,67], it is defined as the ratio of the individual’s wet weight (W, in gram) to the cube of its total length (L_t, in centimeters):

$$\mathrm{F}\mathrm{C}\mathrm{I}=\frac{\mathrm{W}}{{{\mathrm{L}}_{\mathrm{t}}}^3} 100$$

The gonad weight changes according to its reproductive state. To indicate this developmental state, the GSI, which expresses the percentage of gonad weight (W_G, in gram) of the whole-body weight, can be used and was calculated as follows (after [68]):

$$\mathrm{G}\mathrm{S}\mathrm{I}=\frac{{\mathrm{W}}_{\mathrm{G}}}{\mathrm{W}} 100$$

As a measure of the energy status of fish, the HSI is used in terms of the liver weight (W_L) relative to the whole-body weight [69,70]:

$$\mathrm{H}\mathrm{S}\mathrm{I}=\frac{{\mathrm{W}}_{\mathrm{L}}}{\mathrm{W}} 100$$

2.4. RNA Extraction and Hybridization

A subsample of approximately 60 mg of each liver sample was taken under liquid nitrogen to prevent thawing and RNA degradation. Each subsample was homogenized in QIAzol Reagent (Qiagen, Hilden, Germany) and purified following the manufacturer’s instructions. Precipitated RNA was dissolved in 100 µL TE buffer (20 mM Tris/HCl pH 8.0, 0,1 mM EDTA). RNA purity was assessed via UV spectrophotometry (A260/280 NanoDrop, Wilmington, DE, USA) with ratios between 2.02 and 2.13. RNA integrity was further assessed for a randomized subset of samples by capillary electrophoresis (Bioanalyser: Agilent, Waldbronn, Germany). In all tested samples, the RNA integrity was always between 9.4 and 9.8 indicating excellent RNA for the following expression analyses.

Labeling of RNA for microarray analyses was performed with the Low Input Quick Amp Labeling Kit (Agilent) according to the manufacturer’s instructions, with 200 ng total RNA. A 1:16 dilution of a positive control RNA (Agilent RNA Spike-In Kit for one-color arrays) was used to monitor the procedure of sample amplification and microarray workflow. Reference RNA and the RNA samples were labeled with one color, cyanine-3. Labeled and amplified cRNA was purified with the RNeasy kit (Qiagen). The labeling products were quantified using the NanoDrop ND 1000’s microarray measurement protocol.

2.5. Array Design and Hybridization

Prior to the desired 44 K array, a 105 K test-array with 3 probes per sequence of a cDNA-library were designed by Agilent’s eArray-web service to select the best probes for the final microarray. Single-color hybridizations of 4 samples were prepared for 2 runs, 1 representing the gene expression in the liver from sampling round 1 and the other from sampling round 2. Each run comprised the following conditions ([temperature, °C]–[CO₂; ppm]): 0–390, 0–1170, 8–390, 8–1170. The signals from 8 hybridizations (each representing one animal at the respective conditions) were normalized using the R package “limma” (version 3.28.14) and analyzed subsequently. First, cross-correlations of the 8 hybridization signals between the 3 probes per cDNA sequence were calculated, identifying the probe with lowest cross-correlation and excluding it as an outlier. After that, from the remaining list of 2 probes per transcript, the probe showing highest variance across the 8 signals was chosen. Selected probes from pairs with higher than 0.7 Pearson correlations were taken as the best representatives.

The redesigned 44 K microarray was used for 55 samples, 5 samples per combined temperature and PCO₂ treatment. The array included probes for 33.3089 sequences, comprising 3 different sets in total that were optimized on the 105 K test-array: (1) 41.031 probes for 33.179 contigs based on a cDNA library from B. saida (33.179 probes representing 1 contig in the illumina library plus 7852 “questionable probes”) that were based on a cDNA library (total number); (2) 378 probes for 126 sequences of the ion regulatory transcriptome (used in [72]); and (3) 9 probes for 3 pre-selected mitochondrial sequences.

Hybridization of all samples was performed in a standardized manner with respect to volume and yield using the Gene Expression Hybridization Kit (Agilent) according to the manufacturer’s recommendations for the 4 × 44K array format. Dye-labeled samples containing 1.65 µg sample and reference each were hybridized for 17 h at 65 °C. After washing and dye stabilization, slides were directly scanned with an Agilent G2565AA scanner using the AgilentHD_GX_1Color protocol (Agilent Technologies, Santa Clara, CA, USA).

2.6. Statistical Analyses of Gene Expression

The gene expression signals were extracted from the microarray using the R package “limma” (version 3.54.2). Fluorescence signal intensities for each probe on the microarray were extracted, background-subtracted and quantile-normalized between each array in order to obtain similar intensity distributions across the set of arrays [73], and replicate spots were averaged.

The contigs present on the arrays were newly annotated using the high-throughput sequencing module of CLC Genomics Workbench 20 (CLC bio, Aarhus, Denmark) and the latest version of the respective transcriptomes. The sequences were mapped against the Gadus morhua transcriptome (NCBI RefSeq assembly: GCF_902167405.1), resulting in a mapping success of 99% (46413/46857). Afterwards, we filtered the annotations for their e-values with a cut-off of 10⁻²⁰ to ensure correct annotation. This resulted in an annotation success of 80% (37332/46857). Since we specifically designed one probe per contig sequence, duplicate gene annotations were treated as isoforms/genes and kept separate. Later, they were labeled with consecutive numbers according to gene abbreviation (e.g., AGPATA4_1, AGPATA4_2).

We performed a redundancy analysis (RDA) with scaled and centered variables, Hellinger-transformed expression data, and forward selection of the explanatory variables with an alpha threshold of 0.05 using the R package “vegan” (version 2.6-4).

2.6.1. Differentially Expressed Genes

The gene expression of each combined treatment group was compared to the gene expression of the control group of 0 °C and 390 ppm CO₂ to identify differentially expressed genes (DEG) in R using the “limma” package (version 3.54.2, [74,75,76]). “limma” uses linear model features to assess the differential expression between different experimental conditions. Briefly, it models the expression levels of each transcript as a linear combination of experimental factors and covariates. “limma”’s implemented empirical Bayes’ statistics compute moderated t-statistics, moderated F-statistics, and log-odds of differential expression through empirical Bayes’ moderation of the standard errors towards a global value. We filtered the statistical results for adjusted p-values ≤ 0.05, which set the level of significance and indicated genes differentially expressed in comparison to the control.

2.6.2. Hierarchical Clustering

Hierarchical clustering was used to compare DEG patterns across combined treatments. Row-wise clustering of the log2 fold-changes of DEG was based upon a Euclidean distance matrix paired with a Ward agglomerative linkage method. Column-wise hierarchical clustering for treatments was based upon a Euclidean distance matrix paired with a complete linkage method. The analyses were conducted using the packages “gplots” (version 3.1.3) and “stats” (version 4.2.2) package in R [77].

2.6.3. fiNOG Annotation of DEG

To identify the most important transcriptomic changes, we assigned the DEGs to fish-specific orthologous functions (fiNOGs, [78]), which are summarized in categories of orthologous groups (COGs). Therefore, transcripts were assigned with orthologous model sequences from the eggNOG (version 4.0; evolutionary genealogy of genes: Non-supervised Orthologous Groups) database. The microarray sequences were blasted against the fish-specific sequence models using rpstBLASTn (reverse position specific BLAST of nucleotides, NCBI local BLAST tool version 2.2.25+), with an e-value cut-off of 10⁻²⁰.

Based on the fiNOG analysis, we identified the COGs most involved in the transcriptomic response according to the total number of DEGs represented in the COGs. For these DEGs, we used the PANTHER^TM (Protein Analysis THrough Evolutionary Relationships) classification system (version 18.0, [79,80]) to assign them to the gene ontology (GO) terms and create our own GO-slims in order to further specify the transcriptomic response. GO-slims were selected after the list of GO-terms originally assigned to a gene using PANTHER^TM, considering the COG to which it was allocated (for more detail, see Supplementary Table S3). Afterwards, we estimated the percentage contribution of this GO-slim within a gene cluster.

2.7. Statistical Analyses of Physiological Parameters and DEG Patterns

The following tests were accomplished using the R environment version 4.2.2 [77]. We followed the recommendation of evidence-based language [81] for describing our results. The following thresholds were used: p > 0.05 (no evidence of effect), p ≤ 0.05 (moderate evidence), p ≤ 0.01 (strong evidence), and p ≤ 0.001 (very strong evidence of effect).

The data were tested for normal distribution using either the Shapiro–Wilk test or the Kruskal–Wallis test. Homogeneity of variance was tested using Levene’s test. In case of normal distribution and homogeneity, a two-way ANOVA (SGR, body-weight, FCI, HSI, and sexual maturity) was performed to test for temperature- and PCO₂-dependent effects, as well as for interactions between the two parameters. For the datasets with unequal variances, Welch’s ANOVA (GSI, DEG-patterns) was performed. In case of significant effects, a subsequent Student–Newman–Keuls (SGR, body weight, FCI, HSI, and sexual maturity) or Tukey (GSI, DEG-patterns) post hoc test for pairwise comparisons was applied.

3. Results

3.1. Whole-Animal Energetics

3.1.1. Growth during Acclimation

In general, the specific growth rate (SGR) was not significantly affected, either by temperature or CO₂, and no interaction of the two factors was observed. However, the highest PCO₂ increased the variance in SGR at 0, 3, and 8 °C, and may even indicate a threshold level for a decrease in body weight compared to the individual starting weight. The mean SGR accounted for 0.36% body weight gain per day (Figure 1).

3.1.2. Physiological Condition at the End of Acclimation

Fulton’s condition factor (FCI, Figure 2B) was significantly affected by temperature (p ≤ 0.05), with the highest FCI of 0.69 ± 0.007 at 3 °C compared to all other temperature groups (p ≤ 0.001), Correspondingly, the total weight (Figure 2A) was also significantly higher at 3 °C compared to the control (p ≤ 0.01) and highest temperature (p ≤ 0.05). Above 3 °C, weight and FCI decreased with temperature. When comparing between temperature treatments, the weight never fell below the average weight at the control temperature (0 °C), while the FCI continued to decline from 3 °C to both sides of the assessed temperature window. The HSI (Figure 2C) significantly increased with the temperature (p ≤ 0.05) while GSI (Figure 2D) decreased almost linearly with the temperature (p ≤ 0.01) Consequently, the GSI reached its minimum of 1.4 ± 0.5 at the highest temperature.

The ANOVA analysis revealed neither significant influences of PCO₂ on any of the physiological condition parameters nor any interactive effect of temperature and PCO₂ (Figure 2).

3.2. Whole Transcriptome

3.2.1. Gene Annotation

The microarray sequences were annotated according to the Gadus morhua transcriptome (derived from gadMor3.0; GCA_902167395.1), with a mapping success of 80% (37332/46857), after filtering the annotations for the e-value with a cut-off of 10⁻²⁰ to ensure correct annotation.

3.2.2. Redundancy Analysis (RDA)

Within a redundancy analysis (RDA), we identified variables (composed of the individuals’ metadata and water parameters, Supplementary Table S1) as important drivers and effect sizes of transcriptomic changes through a forward selection approach that used permutation of residuals. Our model explained 13% (RDA’s adjusted R²) of the variation in the transcriptomic response across treatments (p ≤ 0.001). The first canonical axis was tested as significant (p ≤ 0.001, RDA2: p = 0.055, RDA3: p = 0.080). The RDA revealed a pronounced temperature effect (p ≤ 0.001). This effect became most obvious by a continuous shift of all ellipsoids of treatments along with the temperature vector in Figure 3, with the 3 °C treatments being in the center of the RDA. Further, a significant influence of the individual fitness of the fish on the transcriptome was detectable, since GSI (p ≤ 0.05) and SGR (p ≤ 0.05) were identified as important effect sizes in the overall transcriptomic response.

3.3. Differentially Expressed Genes (DEGs) and Their Function

3.3.1. DEGs

We analyzed each treatment in comparison to the control condition of 0 °C and 390 ppm CO₂ for differentially expressed genes (DEGs). A log2 fold-change of 0.5 was set as the threshold for all further analyses.

The overall transcriptomic response comprised all significantly up- and downregulated DEGs at each temperature treatment, significantly increasing with temperature (p ≤ 0.001; Figure 4A, Supplementary Table S2). The post hoc test revealed significant differences at 3 °C in comparison to all other temperatures (p ≤ 0.001). At 8 °C, 5071 significant DEGs were found in comparison to the control temperature (p = 0.05). Compared to 3 and 6 °C, this represents a doubling of the genetic response (Supplementary Table S2). Progressive hypercapnia did significantly affect the general transcriptomic patterns in combination with temperature (p ≤ 0.001). At 6 and 8 °C, rising PCO₂ significantly reduced (p ≤ 0.001) the number of DEGs, whereas it triggered the transcriptomic response at 3 °C, causing a four-fold increase in significantly downregulated DEGs and doubling of significantly upregulated DEGs compared to the control PCO₂ (Figure 4B,C, Supplementary Table S2).

In general, the strength of the transcriptomic response (log2 fold-change) within a combined treatment was almost balanced between up- and downregulation (Figure 5, Supplementary Table S2).

After hierarchical clustering, two major clusters became visible (Figure 5, x-axis). The pattern of DEGs assigned the treatments with the highest temperature in combination with the two lower CO₂ concentrations into one cluster that was clearly separated from the remaining treatments. The second cluster could be divided into two sub-clusters according to temperature: One could be described as a warm-temperature cluster, as it contained all three PCO₂ treatments at 6 °C and the high-PCO₂ treatment (1170 ppm) at 8 °C. The second sub-cluster could be described as a low-temperature cluster. Together with the control treatment (0 °C and 390 ppm CO₂), it included all three PCO₂ levels at 3 °C and the highest PCO₂ level at 0 °C.

3.3.2. Fish-Specific Orthologous Functions (fiNOGs)

To identify and categorize the most important transcriptomic changes, we assigned the significant DEGs to fish-specific orthologous functions (fiNOGs), which are divided into categories of orthologous groups (COGs; see also Figure 6 for downregulated COGs and Figure 7 for upregulated COGs). Significant DEGs that were annotated with fiNOGs but belonged to the “poorly characterized” category (general function prediction only: “R”; function unknown: “S”; and no result: “X”) were excluded from further analysis.

Overall, from 4620 significant DEGs, 3227 were annotated with fiNOGs. This corresponds to an annotation success of 70%. In total, 313 of the annotated DEGs were assigned to more than one category. Most DEGs belonged to the category “signal transduction mechanisms” (“T”, 619), followed by “posttranslational modifications, protein turnover, chaperones” (“O”, 566). A total of 1721 of the 3227 significant DEGs were downregulated, of which 142 were assigned to more than one category (Figure 6), and 1506 were upregulated, of which 171 were assigned to more than one category (Figure 7).

In general, most of the DEGs for temperature-induced downregulation were found at 8 °C (Figure 6). An increase in PCO₂ modulated the transcription as well, most obviously at 3 °C, which we identified as the optimum temperature for growth (Figure 1 and Figure 2A,B). At 3 °C and high PCO₂, genes corresponding to signal transduction mechanisms (“T”; 84 DEGs), post-translational modifications, protein turnover, chaperones (“O”; 83 DEGs), and lipid transport and metabolism (“I”; 82 DEGs) were strongly downregulated. However, this PCO₂ effect was lost at 8 °C, probably masked by the overall stronger negative response at this temperature.

A similar picture emerged from analyses of upregulated DEGs, but the PCO₂ effect was not as pronounced as that for downregulation (Figure 7).

Transcriptomic Signatures and Pathway Shifts

The fiNOG analysis revealed that specific molecular mechanisms were substantially influenced by the different combined exposure conditions. Only COG categories with an average finNOG count of more than 90 in at least one treatment were selected, characterizing a strong involvement in the transcriptomic response (see Figure 6 and Figure 7). These comprised (1) lipid transport and metabolism (COG “I”, Figure 8); (2) post-translational modification, protein turnover, and chaperones (COG “O”, Figure 9); (3) amino acid transport and metabolism (COG “E”, Figure 10); and (4) carbohydrate transport and metabolism (COG “G”, Figure 11). Using the PANTHER^TM classification system, we further specified the gene function into GO-slims and estimated the percentile contribution of this GO-slim within a gene cluster to reveal the structures and pathway shifts of the transcriptomic response.

The global transcriptomic response (Figure 5) indicated a strictly temperature-driven differential gene expression, yet when comparing it with the transcriptomic patterns within the four super-ordinated COGs, differentiated pictures emerged: The hierarchical cluster analysis of the DEGs, which belonged to lipid transport and metabolism pathways (Figure 8), resulted in a clear separation of all 6 °C and 8 °C treatment groups from the remaining lower temperatures, which was characterized by a mixed pattern of the respective PCO₂ groups. The hierarchical clustering of DEGs relating to post-translational modifications, protein turnover, and chaperones (Figure 9) provided a very distinct temperature pattern, with only the treatments at 8 °C with medium and control PCO₂ forming a separate cluster, while the other treatments were grouped together. In this large mixed cluster, only 3 °C and high PCO₂ might be considered as a subgroup, emphasizing the impact of hypercapnia at lower temperatures. Interestingly, the separate hierarchical cluster analysis of DEGs related to amino acid transport and metabolism (Figure 10) and carbohydrate transport and metabolism (Figure 11) resulted in another modification of the aforementioned: Again, a clear high temperature cluster comprising all three PCO₂ levels at 8 °C on the one hand and a low-temperature cluster on the other became obvious. But here, the 6 °C/high PCO₂ treatment clustered together with the low temperatures, whereas 3 °C/medium PCO₂ formed a sub-cluster with the remaining 6 °C groups. As observed before, CO₂ seemed to modify the predominant temperature signal for this gene category.

• Lipid Transport and Metabolism

The functional analyses revealed 103 significant DEGs indicative of changes in lipid metabolic processes and transport (Figure 8). The majority of genes represented here were downregulated with increasing temperature. The hierarchical clustering of these genes indicated a temperature-induced decline in lipid biosynthetic processes and signs of impaired energy production by fatty acid catabolism. Only the transcription of genes related to long-chain fatty acid metabolic processes (e.g., diacylglycerol O-acyltransferase 2 (dgat2), mitochondrial medium-chain acyl-CoA ligase (acsf2)), and lipid transport (e.g., fatty acid-binding proteins 3 and 7 (FABP3, FABP7)) was enhanced in treatments above 3 °C (see Figure 8, DEG-cluster 1). Although the number of genes involved in beta-oxidation was not large, their transcription was significantly altered. Interestingly, fish exposed to 3 °C and high PCO₂ showed a strong response in beta-oxidation, as well as when incubated at 8 °C in combination with control and high PCO₂. In particular, the genes ACAA2, encoding the protein 3-ketoacyl-CoA thiolase, which catalyzes the final step of the mitochondrial fatty acid beta-oxidation, and ACOX3, encoding acetyl-CoA oxidase 3, indicated reduced beta-oxidation under increased temperatures. Lipid breakdown and membrane lipid synthesis also contributed to a large percentage of the transcriptomic response (see Figure 8, DEG-cluster 2). Starting at 3 °C and high PCO₂, both processes decreased significantly and to a greater extent with increasing temperature, but above 3 °C, they were independent of PCO₂. In particular, cholesterol synthesis was affected, indicated by a significant decrease in mevalonate phosphate decarboxylase (Mvd), hydroxysteroid 17-beta dehydrogenase 7 (Hsd17b7), and 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR), substantially taking part in the biosynthesis of cholesterol, in most of the treatments above 3 °C [82].

• Post-translational Modification, Protein Turnover, and Chaperones

The transcriptomic responses related to post-translational modification, protein turnover, and chaperones (Figure 9) were more balanced between up- and downregulation compared to lipid metabolic pathways (Figure 8). In total, 161 significant DEGs were related to this functional group. A typical transcriptomic response to temperature exposure close to the thermal maximum of a species would be a pronounced heat shock response in the form of targeted upregulation of heat shock proteins and their additional transcriptomic factors (reviewed by [83]). This did not occur to the extent we expected, as we even found a downregulation of genes encoding for heat shock proteins such as hsp90a, hsp90ab1, HSP90B1, and HSPA5 (see Figure 9 DEG-clusters 3 and 4). Interestingly, this response was most pronounced at 8 °C (in combination with medium PCO₂). In parallel, genes related to hypoxia and oxidative stress were transcribed to a greater extent at this temperature. Besides genes directly related to stress response, genes related to proteolysis, ubiquitination, and protein folding were modulated compared to control conditions. The majority of transcripts related to ubiquitination were downregulated (see Figure 9 DEG-cluster 3), including different proteasome subunits (e.g., proteasome 26S subunit, non-ATPase 12 (PSMD12), proteasome 20S subunit alpha 4 (Psma4), or proteasome 20S subunit beta 3 (psmb3)) and ubiquitin–protein ligases (e.g., E3 ubiquitin–protein ligase rnf168 (rnf168), E3 ubiquitin–protein ligase RNF217 (Rnf217), or ubiquitin-conjugating enzyme E2 Z (ube2z)). Genes involved in protein turnover, e.g., those related to protein folding or proteolysis, were particularly upregulated at the highest incubation temperature.

• Amino Acid Transport and Metabolism

A total of 51 significant DEGs encoding enzymes involved in amino acid transport and metabolism were identified by fiNOG annotation (Figure 10). The transcriptomic response was dominated by an upregulation of catabolic processes and amino acid transport. In particular, catabolic processes of essential amino acids were upregulated at 6 and 8 °C throughout the PCO₂ range (see Figure 10 DEG-clusters 1&2), including genes related to arginine (OAT: ornithine aminotransferase), leucine (AASS: alpha-aminoadipic semialdehyde synthase; Mccc1: methylcrotonyl-CoA carboxylase subunit 1), and tryptophan catabolism (haao: 3-hydroxyanthranilate 3,4-dioxygenase; tdo2a: Tryptophan 2,3-Dioxygenase).

• Carbohydrate Transport and Metabolism

A total of 97 significant DEGs belonged to the functional group of carbohydrate transport and metabolism (Figure 11). Within this group of transcripts, elevated temperature and PCO₂ appeared to induce isoform differentiation, in particular of glucose transporters. Six different glucose transporters (SLC37A4, SLC2A10, SLC2A11, SLC2A2, SLC2A5, and SLC2A9) were differentially expressed, and a total of 28 contigs were found. Most isoforms of the same transcript were expressed in the same directional manner (predominantly upregulated), except for SLC2A2, where some were upregulated while the majority were downregulated. It is also notable that the mammalian acidic chitinases (encoded by CHIA) formed a distinct, highly upregulated gene cluster. Chitinases in the digestive systems of fish hydrolyze the chitin of crustaceans and aquatic insects and play an important role in energy metabolism, but they are also distributed beyond the digestive tract and even play a role in the immune response to pathogens [84,85,86]. An upregulation in liver could therefore indicate some kind of stress-induced response. Except for the 3 °C/normocapnia treatment, all other treatments induced a significant upregulation of the gene. In addition, we found indicators of impaired glycosylation (related genes: ugt3: UDP-glucosyltransferase protein 3, Ugt1a5: UDP-glucuronosyltransferase 1A5, Fut9: fucosyltransferase 9, GCNT7: glucosaminyl (N-acetyl) transferase family member 7), which was apparently induced by temperature early on, as even 3 °C and normocapnia led to different expression of at least some of these genes compared to the control treatment.

4. Discussion

In this study, we correlated the hepatic transcriptomic response of Polar cod with previously measured physiological responses to ocean warming and acidification scenarios [56,57,58,59,60,61]. In this comprehensive project, we conducted our transcriptomic study on a subset of the physiologically tested individuals, providing the invaluable advantage of following the physiological acclimation processes at the transcriptomic level. Furthermore, animals treated at intermediate PCO₂ levels and the respective temperatures were included in the present study to identify possible fine-tuning of increased PCO₂ at the molecular levels, which is not detectable in the former studies at whole-animal and physiological levels. The physiological studies revealed a stronger effect of temperature than PCO₂, which we were also able to track at the transcriptomic level. Although rather modest transcriptomic changes were observed (approximately a mean log2 fold-change of ±1 in all treatments), we were able to identify important pathways for potential acclimation mechanisms in metabolic energy production and maintenance of organismal function.

4.1. General Response to Temperature

Many studies have already postulated a low acclimation capacity of Polar cod (e.g., [46] and references therein, [58]). The modest transcriptomic changes and the lack of a pronounced cellular stress response may indicate that the fish were able to acclimatize to the experimental conditions to a large extent. Nonetheless, a cellular stress response could, of course, have been induced early in the exposure period or during the four month of acclimation and was no longer detectable in the end. Nevertheless, in a similar study on Antarctic eelpout, Pachycara brachycephalum, which was similarly acclimated long-term (−1 °C up to 9 °C), a clear differentiation between successful acclimation and a permanent stress response beyond a tipping point could be identified [87]. Not only were several genes of the cellular stress response [88] expressed at high temperatures, but the overall transcriptomic response was maximal at these extreme temperatures and minimal within the thermal optimum. This led to the conclusion that a new steady state could not be reached beyond the tipping point. Accordingly, even if the Polar cod displayed a cellular stress response in the early phase of the present experiment, this response was not detectable in the long term, so a successful acclimatization with an adapted steady state could be assumed for all applied temperatures (also in combination with hypercapnia).

Furthermore, it also seems that gadoids in general have evolved somewhat different response mechanisms to heat stress than other teleost fish species. In other gadoid species, such as Atlantic cod, Gadus morhua, and haddock, Melanogrammus aeglefinus, the absence of heat shock protein 70 (hsp70) response has been reported [89,90], while individuals have already shown physiological impairments due to warming, such as a reduced growth rate [89]. The lack of a “traditional” heat shock response, which would lead to the identification and removal of damaged proteins, implies that Polar cod has developed a different strategy. This could be covered by a higher protein turnover in general, mediated either by proteolysis or by ubiquitin-mediated degradation in the proteasome, which we found to be upregulated at the highest temperature (Figure 9). Additionally, studies have indicated that a profound heat shock response is absent in some cold-adapted species, such as the emerald rockcod, Trematomus bernacchii, or the black cod, Notothenia angustata [91,92,93]. This may be because they never encounter heat waves and have lost it evolutionarily (to save energy), or, as Buckley et al. (2009) [91] suggested, because the cells of Antarctic fish and other cold-adapted species may already be operating at their maximum capacity to produce heat shock proteins, since cold shock reactions may play a role in protein folding similar to that of heat shock. High levels of constitutive hsp expression may therefore attenuate an early onset of heat shock response. Short-term exposures at critical temperature for Polar cod would also be necessary to determine the onset of classic heat stress reaction. Nonetheless, the present transcriptomic responses in conjunction with the physiological parameters clearly indicate no severe or chronic warming stress for Polar cod even after four months of exposure to 8 °C, which is far beyond their natural habitat temperatures.

4.2. General Response to PCO₂

Our study brought small, but significant, PCO₂-induced impairments to light, which were also reflected in the whole-animal studies. A temperature-independent reduction in behavioral laterality was found, suggesting reduced fitness of the Polar cod under future CO₂ conditions [60]. This was corroborated by an impaired swimming performance under elevated PCO₂, also indicating insufficient oxygen turnover, as the oxygen demand could not be met despite the maximum metabolic rate being increased [57]. We also found CO₂ effects, even though most of them were interactive with temperature and thus difficult to specify on the process levels. This CO₂ effect was most obvious at 3 °C; where the progressive rise in PCO₂ significantly increased the number of differentially expressed genes (Supplementary Table S2, Figure 9, Figure 10 and Figure 11). Interestingly, the combined stressors showed no influence at the edges of the temperature spectrum, but did so in the range of the species’ optimum temperature. In this range, the temperature effect appeared to be small enough that a PCO₂ effect became detectable. Even though the oceans’ CO₂ levels will rise in parallel with increasing temperatures, these cold but high-CO₂ conditions will occur in the natural habitat of the Polar cod if global warming continues to progress. Polar cod is known to inhabit fjord systems, where the loss of sea ice already increased summer stratification in deep fjords [94,95]. An extension of ice-free winters for more than one season enables the formation of low-oxygen and CO₂-rich conditions in the cold deep-water layers, as these can no longer be supplied by sinking oxygen-rich surface water created during ice formation [94,96,97], a process referred to as thermohaline circulation [98]. It can therefore be assumed that Polar cod will experience considerable impairment due to an increased PCO₂ level at relatively cold temperatures.

4.3. Whole-Animal Effects and Metabolic Adjustments

Compared to the results of Kunz et al. (2016) [56], who investigated the growth performance and feed conversion of Polar cod in our acclimatization experiment (n = 12 per combined treatment), the subgroup of fish we used for our transcriptome study (n = 5 per combined treatment) showed a slightly lower temperature for optimal growth (SGR maximum: 3 °C, and 4.7 °C in a fitted model for maximum FCI). However, if we also take the HSI and GSI into consideration, it becomes clear that the optimal temperature for general development also appears to be somewhat slightly higher than the temperature for optimal growth. The HSI reached its maximum at the highest temperature, while GSI decreased with temperature and plateaued between 3 and 6 °C, and temperatures above 6 °C impaired gonad maturation. Additionally, our subgroup consisted of only male individuals. Studies have shown that exposure to high temperatures during development can lead to reduced reproductive output and poorer-quality offspring in males, while females may experience enhanced reproduction and produce higher-quality offspring under similar conditions [99]. Another sex-specific response to temperature can occur in some fish species during sex determination, highlighting the different sensitivities of the sexes to changing environmental conditions. Potentially through epigenetic mechanisms in response to temperature, genotype–phenotype mismatches during development could lead to long-term impacts on population dynamics such as masculinization or feminization [100,101,102]. Furthermore, altered water temperature regimes, particularly from a physiological optimum to higher temperatures, can influence gonadal development in fish, impacting cell division, nutrient accumulation, or GSI differently in males and females [103,104]. Indeed, the physiological studies on the entirety of the incubated fish (n = 12, including males and females [56]) revealed differences in gonadal maturation between the sexes. The variability of GSI tended to decrease in male Polar cod with increasing temperatures, while females generally showed a low GSI over the entire temperature range [105]. These findings collectively suggest a differential susceptibility of male and female gonad production to elevated temperatures in fish. For our small, strictly male subgroup (n = 15 per temperature (n = 10 for 0 °C)), we would therefore assume an optimal temperature for growth and development between 3 and 6 °C, which also fits to the species’ preferred temperature for activity and metabolic capacity reported by Schurmann et al. (1994) [106] and Leo et al. (2017) [58].

Different studies have already examined the acclimatization potential of Polar cod to temperature, and have indicated maximum growth at temperatures between 5–6.5 °C [56,107,108] and a general preferred temperature for activity and metabolic capacity between 3–6 °C [58,106]. From 6 °C onwards, we detected important metabolic shifts in the hepatic transcriptome. Our findings indicated a warming-induced decline in lipid biosynthetic processes and signs of reduced energy production by fatty acid catabolism, whereas carbohydrate metabolism could be largely maintained. The aerobic respiration quotient of carbohydrate, fatty acid, and protein metabolism indicates that carbohydrates have the lowest oxygen demand and also function at limited oxygen supply capacity, as is postulated at warm temperatures [24]. Since carbohydrate metabolism supports anaerobic metabolism through glycolysis, it is possible that this metabolic pathway is favored at warmer temperatures. This strategy has also been observed in another cold-adapted fish species, such as the Antarctic eelpout Pachycara brachycephalum [87,109], and the present results for Polar cod might indicate hypoxemia at warmer temperatures and point to the same adaptive strategy as in highly cold-adapted Antarctic fish.

4.4. Carbohydrate Metabolism

Although earlier studies have revealed that Polar cod has rather little anaerobic metabolic capacity ([57,62] under review), we found signs of a decrease in cellular oxygen saturation with increasing temperature at the transcriptomic level, as well as a preparation of the animals for a shift to metabolic pathways that can withstand functional hypoxemia, which is more likely to occur with warmth.

Genes related to hypoxia were upregulated in warmer treatments such as the endoplasmic reticulum oxidoreductin-like protein 1α (ero1a). Decreasing cellular oxygen tension stimulates upregulation of the expression of endoplasmic oxidoreductase L α [110]. An oxygen-induced regulation of ero1a is supposed to maintain the transfer rate of oxidizing equivalents to the protein disulfide isomerase, which is essential for the formation of disulfide bonds in the endoplasmic reticulum [110]. A number of endoplasmic proteins are activated by a malfunctioning energy production system, induced, for example, by severe cellular hypoxia or glucose deprivation [111]. A sign of glucose deprivation might be the extensively upregulated transcription of glucose transporters and the large number of isoforms that are present. Some are referred to as bidirectional glucose transporters, some strictly as importers. Since the liver is an important key organ for supplying the whole organism and the different tissues with energy resources, the data indicate increased glucose turnover under warming. Environmental factors like hypoxia and hypercapnia can influence the expression of glucose transporter genes [112,113]. For instance, in sea bass (Dicentrarchus labrax), hypoxia lead to an increase in hepatic GLUT2 mRNA levels [112], while in spiny chromis (Acanthochromis polyacanthus), different glucose transporters such as gtr1 (encoded by SLC2A1) show transgenerational upregulated expression under elevated PCO₂ [113]. Additionally, studies on zebrafish (Danio rerio) and coho salmon (Oncorhynchus kisutch) indicate that elevated growth hormone levels could modulate various glucose transporter expressions, affecting energy uptake [114,115]. Therefore, it is plausible that increased temperature, as an environmental driver, could potentially affect the expression of genes related to glucose transporters in fish similar to the effects observed under hypoxic or hypercapnic conditions and altered hormone levels. Hence, this upregulation may be considered a stress-related response mechanism.

4.5. Mitochondrial Functioning

The physiological parameters indicate a disrupted balance between energy production and demand, leading to a shift in energy allocation. At 8 °C, ATP production was impaired and an increased proton leak and decreased complex IV capacity was observed [58], while the metabolic rate was increased [57]. We were able to trace signs of mitochondrial dysfunction at the transcriptomic level in the liver: The mitochondrial protein peroxiredoxin 3 (PRDX3), a protein with antioxidant function, was upregulated at the highest temperature. Peroxiredoxins (Prxs) act as defensive molecules to reduce the excessive amounts of reactive oxygen or nitrogen species (ROS/RNS) and regulate their mediatory role [116]. Peroxiredoxins are also found in other fish, such as rock bass (Oplegnathus fasciatus) and Antarctic emerald rockcod (Trematomus bernacchii), in response to environmental warming and/or hypoxia [117,118]. Another protein repair system activated by the presence of ROS is the methionine sulfoxide reductases system (methionine sulfoxide reductases A or B). We found two representatives of this enzymatic family to be upregulated with temperature: MSRB2 and 3. These enzymes catalyze the reduction of methionine sulfoxide, which is produced in the oxidative process of methionine side chains in the presence of ROS (reviewed by [119]). Different environmental drivers can lead to oxidative stress, such as thermal stress [120,121]. However, small amounts of ROS are also synthesized naturally under normal conditions. Palma et al. (2024) [122] found that over 98% of the consumed oxygen is used for ATP production, while only 1–2% is converted to ROS. For a long time, ROS were associated with mitochondrial dysfunction and disease (as reviewed by [122]). Although recent studies have proven ROS production to be tightly regulated and coupled to an activation of adaptive processes, and to play important roles in cell signaling and gene expression (reviewed by [122]), we would agree with the common interpretation and attribute the presence of ROS in our organism to mitochondrial dysfunction at high temperatures.

4.6. Lipid Metabolism

Among the four molecular mechanisms of energy production that we have studied in detail, lipid metabolism appears to be the most sensitive to changes in environmental conditions. In particular, lipid degradation and membrane lipid synthesis decreased significantly from 3 °C and high PCO₂ with increasing temperature.

The homeoviscous adaptation of membrane fluidity [123] is a well-established acclimation concept to regulate the viscosity of membrane lipids in response to changing temperature conditions. The relationship between membrane fluidity and temperature is generally inverse [124,125]. Organisms usually increase their membranes’ fluidity when temperatures are decreasing and vice versa. In general, cholesterol increases membrane stability at higher temperatures and is usually found in higher concentrations at higher temperatures [126,127]. Interestingly, we detected signs of increased membrane fluidity and specifically a decrease in stabilizing mechanisms such as cholesterol with increasing temperatures. The parallel steep increase in long-chain fatty acid metabolism at higher temperatures might indicate that Polar cod increased their membrane fluidity by implementing more long-chain fatty acids in their membranes in the warm treatments. Polyunsaturated long-chain fatty acids are especially involved in a network of pathways in fish; they are known to influence enzyme activities and gene expression related to lipid metabolism pathways, to regulate membrane fluidity, as well as to play a role in the immune and stress response [128,129,130]. Unfortunately, the transcriptomic data do not provide any indication of whether the saturation state of the long-chain fatty acids was affected. Besides these findings on hepatic membrane modulation, Leo et al. (2020) [59] did not detect signs of a modulation of cardiac cell membrane composition in Polar cod due to temperature or PCO₂, indicating tissue specific responses to changing environmental conditions.

Interestingly, HSI was the only fitness parameter that increased with temperature, reaching its maximum at 8 °C, although the energy demand increased with temperature and the feed intake decreased significantly [56]. Increasing HSI might indicate that liver homeostasis (and lipid storage [131] was prioritized with progressive warming even though feed conversion was reduced [56]. At the same time, the GSI decreased steadily, indicating that less energy was available for gonad maturation. For Atlantic cod, indications of a positive interaction of increasing temperature and PCO₂ on the liver weight were found [56]. Kunz et al. (2016) [56] suggested that the trend of HSI to increase with temperature might be attributed to a decreasing ability to exploit liver energy reserves in warm conditions. Although the number of genes involved in beta-oxidation was not large, their transcription was significantly reduced upon warming, which would support the assumption of Kunz et al. (2016) [56]. Likewise, fish exposed to 3 °C and high PCO₂ showed a strong negative response in beta-oxidation, as well as when incubated at 8 °C in combination with high PCO₂.

Additionally, the increase in liver weight may also be supported by an increased metabolism of long-chain fatty acids. In many organisms, triacylglycerols are the main stores of metabolic energy (cf. [132], and reviewed by [133]), including fish [134]. Triacylglycerols serve as reservoirs of essential and non-essential fatty acids. The final step of their formation is catalyzed by diacylglycerol O-acyltransferase 2, forming an ester bond between 1,2-diacylglycerol and a long-chain fatty acyl-CoA [135,136,137]. The formation of long-chain fatty acyl-CoA itself is connected to acsf2, belonging to the family of acyl-CoA synthases, which catalyze the initial reaction in fatty acid metabolism [138]. It is not very common that the biosynthesis of long-chain fatty acids is enhanced by an elevation of temperature (e.g., [139,140]); however, it is conceivable that the reduced beta-oxidation also resulted in an increased biosynthesis of long-chain fatty acids. The reason for such an accumulation, apart from the impending energy limitations, must be investigated in future studies. It is possible that these lipid stores will serve as the first source of energy when the conditions for the organism become more favorable again.

4.7. Protein Degradation and Mobilization of Essential Amino Acids

An important finding of this study is the suppressing effect of the highest temperature on the protein degradation. Most affected was the ATP-dependent ubiquitin–proteasome pathway at 8 °C. Different proteasome subunits (e.g., proteasome 26S subunit, non-ATPase 12 (PSMD12), proteasome 20S subunit alpha 4 (Psma4), or proteasome 20S subunit beta 3 (psmb3)) and ubiquitin–protein ligases (e.g., E3 ubiquitin–protein ligase rnf168 (rnf168), E3 ubiquitin–protein ligase RNF217 (Rnf217), or ubiquitin-conjugating enzyme E2 Z (ube2z)) were downregulated. Reduced protein degradation might serve as a compensatory mechanism to conserve energy. This was also reported for other fish species, such as Atlantic salmon (Salmo salar) [141] and rainbow trout (Oncorhynchus mykiss) [142,143] incubated in conditions of either chronic heat stress or temperatures close to the upper thermal limit. At the same time, we also observed increased proteolysis at the highest temperature, highlighting a shift in the energy distribution due to the increased temperature, as proteolysis enables the release of amino acids to synthesize new proteins needed for survival or support energy production as intermediates of the TCA (cf. [144,145]).

Moreover, we found indications of increased degradation of essential amino acids (Figure 10). Various studies have already shown that increased water temperatures lead to changes in amino acid metabolism, including increased degradation of some essential amino acids, for example, in zebrafish (Danio rerio) (tryptophan, methionine) [146] and in the Antarctic marbled rockcod (Notothenia rossii) (leucine, methionine) [147]. In addition, warming and acidification negatively affect the levels of essential amino acids in marine primary producers, potentially impacting higher trophic levels that rely on these essential biomolecules [148]. These results suggest that the degradation of essential amino acids at the upper thermal limit of Polar cod could be seen as a kind of last-resort energy supply, but one that could affect their overall metabolism and health.

5. Conclusions

The hepatic transcriptomic response indicated strong adjustments of pathways related to the central energy metabolism with progressing temperature. Together with whole-animal and physiological parameters, e.g., behavioral impairments, increased energy turnover, increased mitochondrial proton leak, and decreased growth, these data indicate a tipping point of around 8 °C for the long-term performance of the Svalbard population of Polar cod. Temperature appeared to be the more critical parameter compared to ocean acidification. We found only modest transcriptomic changes, revealing a greater acclimatization potential for Polar cod than other studies have previously suggested. This is also supported by the lack of a classical heat shock response, no thermally induced excessive protein degradation, a maximum FCI at 4.7 °C (from the modulated curve), and a clear increase in HSI up to the maximum temperature applied. Nonetheless, a sharp decline in gonadal maturation and GSI suggests that total maintenance costs exceed energy supply to such an extent that they can only be met at the expense of an overall reduction in animal fitness [cf. [24]]. Transcriptomic investigations of further tissues obtained within this study would help to identify other important mechanisms involved.

We identified the following essential acclimatization strategies of Polar cod with progressive warming and/or ocean acidification: (1) As temperature rises, most parts of the lipid metabolism, especially beta-oxidation, were significantly reduced, while (2) the organisms were able to largely maintain their carbohydrate metabolism. (3) The breakdown of essential amino acids appeared to be a kind of emergency-energy-supply strategy at the species’ upper thermal limit, and (4) the downregulation of ATP-dependent ubiquitination possibly saved energy for other processes. These processes may include (5) the accumulation of storage lipids, as indicated by a steep increase in long-chain fatty acid biosynthesis and a parallel increase in HSI with progressive temperature.

Our results, therefore, suggest that Polar cod’s survival strategy for higher-than-optimal temperatures seems to be rather passive. The organism seems to prepare itself to wait for the situation to improve; to store and spare energy in case the situation deteriorates further; and to further limit functions such as reproduction, growth, swimming ability [57], and mobility [60]. However, the question remains as to whether these adaptation mechanisms can be realized in their natural habitat, especially if food supply is limited (in contrast to the present experiment) and at the same time predation pressure on Polar cod will increase as boreal species migrate into the warming Arctic ecosystems.