The transcriptomic responses of Atlantic salmon (Salmo salar) to high temperature stress alone, and in combination with moderate hypoxia

Background Increases in ocean temperatures and in the frequency and severity of hypoxic events are expected with climate change, and may become a challenge for cultured Atlantic salmon and negatively affect their growth, immunology and welfare. Thus, we examined how an incremental temperature increase alone (Warm & Normoxic-WN: 12 → 20 °C; 1 °C week− 1), and in combination with moderate hypoxia (Warm & Hypoxic-WH: ~ 70% air saturation), impacted the salmon’s hepatic transcriptome expr\ession compared to control fish (CT: 12 °C, normoxic) using 44 K microarrays and qPCR. Results Overall, we identified 2894 differentially expressed probes (DEPs, FDR < 5%), that included 1111 shared DEPs, while 789 and 994 DEPs were specific to WN and WH fish, respectively. Pathway analysis indicated that the cellular mechanisms affected by the two experimental conditions were quite similar, with up-regulated genes functionally associated with the heat shock response, ER-stress, apoptosis and immune defence, while genes connected with general metabolic processes, proteolysis and oxidation-reduction were largely suppressed. The qPCR assessment of 41 microarray-identified genes validated that the heat shock response (hsp90aa1, serpinh1), apoptosis (casp8, jund, jak2) and immune responses (apod, c1ql2, epx) were up-regulated in WN and WH fish, while oxidative stress and hypoxia sensitive genes were down-regulated (cirbp, cyp1a1, egln2, gstt1, hif1α, prdx6, rraga, ucp2). However, the additional challenge of hypoxia resulted in more pronounced effects on heat shock and immune-related processes, including a stronger influence on the expression of 14 immune-related genes. Finally, robust correlations between the transcription of 19 genes and several phenotypic traits in WH fish suggest that changes in gene expression were related to impaired physiological and growth performance. Conclusion Increasing temperature to 20 °C alone, and in combination with hypoxia, resulted in the differential expression of genes involved in similar pathways in Atlantic salmon. However, the expression responses of heat shock and immune-relevant genes in fish exposed to 20 °C and hypoxia were more affected, and strongly related to phenotypic characteristics (e.g., growth). This study provides valuable information on how these two environmental challenges affect the expression of stress-, metabolic- and immune-related genes and pathways, and identifies potential biomarker genes for improving our understanding of fish health and welfare. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07464-x.


Background
Temperature and oxygen are key environmental factors that influence the physiology, metabolism and survival of marine organisms, including fish [1][2][3][4][5][6]. Aquatic environments are characterized by short-(i.e., heat waves) and long-term (i.e., seasonal) fluctuations in water temperatures, which may become a further challenge for marine fish species with global warming [7][8][9]. For example, global ocean temperatures are projected to increase by an additional 1-3°C by the end of the century [9], and this will be associated with more widespread and severe periods of hypoxia (low dissolved oxygen [DO]) in coastal regions [10,11].
Thermal stress responses in fish have been widely investigated, and involve the expression of evolutionarily conserved genes [12]. A number of studies have demonstrated that acute and/or long-term exposure to high temperatures results in extensive changes in gene transcription in different salmonid tissues [13][14][15][16][17][18][19][20]. These molecular responses include alterations in the expression of genes related to the heat shock response, protein folding and repair, stress-induced cell death/apoptosis, signal transduction, oxidative stress, the inflammatory response and a diversity of metabolic processes [13][14][15][16][17][20][21][22][23]. Similarly, hypoxia has profound effects on a broad range of biochemical, physiological and behavioural processes in fishes, and has deleterious impacts on growth and reproduction that eventually influence health, welfare and survival [2][3][4]24]. Under hypoxic conditions, fish suppress energy-requiring processes like protein synthesis, aerobic metabolism and mitochondrial energy production [25][26][27][28]. On the contrary, hypoxia stimulates anaerobic ATP production, lysosomal lipid trafficking/ degradation, the antioxidant system, the cellular heat shock response and immune-related pathways [25][26][27][29][30][31]. In previous experiments fish have generally been exposed to an acute temperature increase or constant high temperatures [13-15, 17, 19, 22, 32], rather than the long-term incremental rise in temperature that is seen at aquaculture sites in coastal regions in temperate zones [33,34]. Given the predicted increase in these two environmental stressors with global climate change [7,[9][10][11], it is of great importance that we assess their combined effects in experiments that simulate realworld conditions [35].
The Atlantic salmon is the most important commercially farmed salmonid species in the world. Juvenile and adult salmon reared in sea-cages in the North Atlantic are facing surface water temperatures up to 18-20°C for extended periods in the summer [33,34], while in Tasmania water temperatures inside the cages have already reached~23°C [36]. Yet, Atlantic salmon have an optimal growth performance at water temperatures between 10 and 14°C [37,38]. In addition, water oxygen levels within the cages fluctuate substantially due to temperature, fish density, feeding and low water exchange [39][40][41], and hypoxic events (~60-70% air saturation) often occur in late summer [33,34,36]. These suboptimal conditions may negatively affect the salmon's physiological and growth performance [1,42], and recently led to mass mortalities at cage-sites in Newfoundland, Canada [33]. Consequently, these conditions are raising concerns worldwide about the profitability of the industry and salmon welfare and survival [1,36]. Nevertheless, we have limited knowledge about the capacity of Atlantic salmon to tolerate heat stress in combination with hypoxia, and whether these stressors interact synergistically or antagonistically, or impose additive effects [4,35].
In this study, we explored the hepatic transcriptional response of post-smolt salmon exposed to an incremental increase in temperature (12 → 20°C at 1°C week − 1 ) and normoxia (∼100% air saturation) (Warm & Normoxic-WN) or in combination with moderate hypoxia (~70% air sat.) (Warm & Hypoxic-WH), as compared to control conditions (12°C, normoxia) (Control-CT) (Fig. 1). The WH condition simulates the environmental challenges that farmed Atlantic salmon can experience during the late summer/fall in sea-cages in the North Atlantic [33,34,41]. The liver was chosen to study due to its roles in a number of biological processes including the stress response, nutrient metabolism and immunity [24,43], and because it has previously been shown to be an excellent organ for characterizing temperature and hypoxia stress responses in this species [14,16]. An Agilent 44 K salmonid oligonucleotide microarray platform [44] was employed to initially assess hepatic transcriptome changes, and to elucidate the processes and mechanisms involved in the liver's response once temperature reached 20°C. Further, we performed real-time quantitative polymerase chain reaction (qPCR, Fluidigm Biomark™) on 41 target genes to: i) validate the microarray results; ii) examine the salmon's molecular stress and immune responses to these environmental stressors; iii) correlate gene expression responses to physiological and growth parameters; and iv) identify genes that have potential as biomarkers for improving our understanding of Atlantic salmon health and welfare under sea-cage conditions predicted to accompany climate change, and for potential incorporation into broodstock selection programs.

Significance analysis of microarrays (SAM)
Based on 17,072 detected microarray probes, the pairwise comparisons computed within SAM recognized 1900 differentially expressed probes (DEPs) for WN challenged fish and 2105 DEPs for WH treated fish as compared to CT fish (Fig. 2a). The complete annotation of SAM-identified DEPs for the WN and WH treatment groups is listed in Additional file 1. When comparing the 'WH vs CT' and 'WN vs CT' DEP lists, we found that 1111 DEPs (38%) were overlapping (i.e., were in common), that 789 DEPs (28%) were WN-specific, and that 994 DEPs (34%) were WH-specific (Fig. 2a). In summary, we identified a unique set of 2894 DEPs when considering WN and WH conditions together, out of which 1267 DEPs were up-regulated (WH = 732, WN = 135; shared for WH and WN = 400) and 1627 DEPs were down-regulated (WH = 262, WN = 654; shared for WH and WN = 711) ( Fig. 2a; Additional file 1). A hierarchically clustered heat map of the 2894 DEPs displayed a robust cluster containing all CT samples that grouped separately from the cluster containing both treatment groups, with distinctive opposing patterns of up-and down-regulation (Fig. 2b). According to Ward's cluster algorithm, the profiles of WN and WH fish created a uniform cluster with a similar magnitude of expression (Fig. 2b).  The Principal Component Analysis (PCA) based on 17,072 detected microarray features displayed a similar global transcription profile for fish exposed to WN or WH conditions in comparison to CT fish that was separated along the first Principal Component (PC-1: p = 2e-4, explaining 16.6% of the variance; CT vs WN p = 2e-4, CT vs WH p = 0.001, WN vs WH p = 0.489; Fig. 2c and Table 1). Similarly, we found a comparable global transcript expression response for the shared 1111 DEPs, as shown by equal cluster distribution for the WN and WH treatment groups as compared to the CT group (PC-1: p = 1e-4, explaining 61.5% of the variance; CT vs WN p < 0.0001, CT vs WH p = 2e-4, WN vs WH p = 0.841; Fig. 2d and Table 1). However, the 789 WN-specific DEPs displayed stronger differential expression in WN fish as compared to WH fish, as indicated by a clear cluster separation between all three treatment groups (PC-1: p = 2e-4, explaining 46.6% of the variance; CT vs WN p = 6e-4, CT vs WH p = 2e-4, WN vs WH p = 0.012; Fig. 2e and Table 1); and the 994 WH-specific DEPs displayed stronger differential expression in WH fish, as evidenced by an extreme shift away from the CT and WN groups (PC-1: p = 1e-4, explaining 42.1% of the variance; CT vs WN p = 0.002, CT vs WH p < 0.0001, WN vs WH p = 0.007; Fig. 2f and Table 1).

GO/pathway term network analysis
Salmon exposed to either the WN or WH treatments had similar patterns for enriched 'Gene Ontology' (GO) terms, 'Kyoto Encyclopedia of Genes and Genomes' (KEGG) and 'Reactome' pathways (hereinafter referred to as 'GO/pathway terms') for up-and down-regulated differentially expressed genes (DEGs) (Figs. 3 and 4; Additional files 2-6). For both treatment groups, the upregulated DEGs were functionally associated with significantly enriched GO/pathway terms (p = 1e-2 to p = 1e-4) connected to Heat Shock Response ( # 1), Cellular Stress ( # 2), Oxidative Stress ( # 3), Apoptosis ( # 4), Immune Response ( # 5), Protein Processing & Localization ( # 6) and Transcription ( # 7) (Figs. 3a and 4a; Additional files 2a-b, 3 and 4). On the contrary, the down-regulated DEGs for WN and WH fish were associated, and with a higher significance, to enriched GO/pathway terms (p = 1e-5 to p = 1e-27) related to Oxidative Stress ( # 3), Proteolysis ( # 8), Catabolic Processes ( # 9) and Cellular Metabolic Processes ( # 10), and formed complex and interconnected functional clusters (Figs. 3b and 4b; Additional files 2cd, 5 and 6). Below we list the shared and dissimilar most significantly enriched GO/pathway terms (with identifiers) for WN and WH fish as visualized in the functional ordered networks (Figs. 3 and 4). However, some terms were abbreviated for simplification purposes, and the complete lists of all terms are given in Additional files 3-6.
GO/pathway term networks associated with up-regulated DEGs Heat shock response ( # 1) In the functionally grouped network analyses, WN and WH fish had up-regulated DEGs associated with enriched GO/pathway terms related to Heat Shock Response ( # 1) (Figs. 3a and 4a). WN fish had one large cluster (~8 terms) with the following enriched GO/pathway terms: 'STIP1(HOP) binds HSP90 and HSP70:HSP40:nascent protein' (R-HSA:3371503), 'ATP hydrolysis by HSP70' (R-HSA:3371422); and 'response to topologically incorrect protein' (GO:0035966) ( Fig.  3a; Additional file 3). On the contrary, WH fish had two large and highly interconnected clusters (~30 terms) associated with the chaperone-mediated heat shock response that included the following enriched Reactome/KEGG pathways: 'ATP binding to HSP90 triggers conformation change' (R-HSA:5618107), 'dissociation of cytosolic HSF1:HSP90 complex' (R-HSA:5324632), 'STIP1(HOP) binds HSP90 and HSP70:HSP40:nascent protein' (R-HSA:3371503), 'protein processing in endoplasmic reticulum' (KEGG: 04141), and 'unfolded protein response (UPR)' (R-(See figure on previous page.) Fig. 2 Hepatic transcriptome responses of Atlantic salmon exposed high temperature stress alone, or combined with hypoxia. a Results of Significance Analysis of Microarrays (SAM) with a False Discovery Rate (FDR) of < 5% (see Additional file 1). The top Venn diagram illustrates the total number of differentially expressed probes (DEPs) in salmon exposed to Warm & Normoxic (WN: 20°C, 100% air saturation) or Warm & Hypoxic (WH: 20°C,~70% air sat.) conditions as compared to the Control group (CT: 12°C, 100% air sat.) (n = 6 per treatment, N = 18 total). The bottom Venn diagram shows the corresponding number of DEPs that were up-or down-regulated. The overlapping areas represent shared DEPs between the WN and WH treatment groups. b Hierarchically clustered heat map based on 2894 DEPs (FDR < 5%) using Ward's minimum variance method. Heatmap and dendrograms illustrate the clustering structure of 2894 SAM-identified DEPs for fish subjected to control conditions (CT), high temperature and normoxia (WN), or high temperature and hypoxia (WH). The integrated colour code shows up-regulated DEPs in red and down-regulated DEPs in blue, and represents normalized log 2 ratios. c Principal Component Analysis (PCA) based on all detected 17,072 microarray probes of fish exposed to CT, WN and WH conditions; d PCA based on a common set of 1111 DEPs shared between the WN and WH treatment groups; e PCA based on 789 DEPs in the WN treatment group; f PCA based on 994 DEPs in the WH treatment group. Each PCA plot includes the three different treatment groups (CT, WN and WH), and is based on normalized log 2 ratios. Ellipses denote the dispersion of variance with 95% confidence intervals around the the center of the distribution for each treatment group. The variance explained by the main Principal Components (PC-1 and PC-2) are indicated in percentage values next to the axes Cellular stress ( # 2) and oxidative stress ( # 3) Fish exposed to WN conditions had larger clusters of up-regulated DEGs associated with Cellular Stress ( # 2) and Oxidative Stress ( # 3), and contained the following significantly enriched GO-terms: 'response to stress' (GO:0006950), 'response to oxidative stress' (GO: 0006979), 'response to toxic substance' (GO:0009636) and 'regulation of response to stress' (GO:0080134)

Apoptosis ( # 4)
For WN fish, we found up-regulated DEGs associated with six large clusters containing~92 highly interconnected enriched GO/pathway terms related to Apoptosis ( # 4), such as 'intrinsic apoptotic signaling pathway in response to nitrosative stress' ( Immune response ( # 5) Pronounced differences between the WH and WN groups were found for pathways related to immune response (  (See figure on previous page.) Fig. 3 Functionally grouped gene network analysis for Atlantic salmon exposed to Warm & Normoxic conditions. Network based on a 377 upregulated and b 798 down-regulated differentially expressed genes (DEGs) in salmon exposed to Warm & Normoxic (WN: 20°C, 100% air saturation) conditions as compared to the Control group (CT: 12°C, 100% air sat.) (n = 6 per treatment, N = 12 total). GO-terms and pathways were obtained through a functional enrichment analysis using the ClueGO plugin in Cytoscape (v3.5.1) [45]. Each node represents a significantly enriched Gene Ontology (GO), KEGG or Reactome pathway (hypergeometric test p < 0.05 with Benjamini-Hochberg correction). The node colour regime visualizes functional groups and processes that share similar genes. The most significant enriched terms of each functional group are illustrated as a summary label. The node circle diameter corresponds to the significance of the enriched pathway (i.e., larger circles correspond to higher significance). The gating and numbering represent the classification of ten functional themes: Heat Shock Response ( # 1), Cellular Stress ( # 2), Oxidative Stress ( # 3), Apoptosis ( # 4), Immune Response ( # 5), Protein Processing & Localization ( # 6), Transcription ( # 7), Proteolysis ( # 8), Catabolic Processes ( # 9) and Cellular Metabolic Processes ( # 10). GO/pathway term annotations were obtained using the GO database for biological process, cellular component, molecular function and immune processes, and the KEGG and Reactome pathway databases. Abbreviations of GO/pathway terms were used for simplification, and the complete list of terms is represented in Additional files 3 and 5

qPCR validation of microarray approach
We chose a subset of 41 microarray-identified genes of interest (GOIs) to confirm the microarray results using a gene-targeted qPCR approach ( Table 2). A significant positive correlation (R = 0.87; p < 1.7e-13) between the mean log 2 fold change (FC) values of the 41 GOIs measured via the 44 K microarray and qPCR methods indicates that the microarray results were validated with great confidence (Additional file 7). In a direct comparison between mean FC values of the microarray and qPCR approaches, 34 genes (83%) showed similar transcription patterns and FC values (Additional file 8). Minor differences were found for the genes camp-a, cat and cldn3, while igfbp2b1, irf2, nckap1l and tapbp had a difference in the direction of expression change (Additional file 8). Out of the 41 GOIs, we found 23 with a significant 'treatment' effect in the linear mixed-effect model (LME), while the other 18 genes were differentially expressed in the SAM analysis but not validated using qPCR (Table 3). This discrepancy may have been caused by: i) using two different technologies; ii) a different sample size; iii) two different statistical approaches (SAM vs LME); or iv) paralog crosshybridization on the 44 K array. The relative expression values [i.e., relative quantity (RQ)] for a selection of GOIs are shown in Fig. 5.

Gene-by-gene analysis
Genes related to heat shock response ( # 1) Expression of the heat shock response-related genes ser-pinh1 (p = 0.010; Fig. 5a) and hsp90aa1 (p = 0.006; Fig.  5b, Table 3) was higher in fish challenged with the WN and WH conditions in comparison to the CT group. On the contrary, while the gene hsp70 was slightly elevated in both treatment groups, this effect was not significant (p = 0.181; Fig. 5c, Table 3).
Genes related to cellular stress ( # 2), oxidative stress ( # 3), transcription ( # 7) and cellular metabolic processes ( # 10) The stress-related gene hcn1 was significantly upregulated in WN and WH fish as compared to CT fish (p = 0.006; Fig. 5d, Table 3). In contrast, seven other stress-related genes were significantly down-regulated in both treatment groups as compared to the CT group: cirbp (p = 0.002; Fig. 5e), calm (p = 0.002; Fig. 5f), cyp1a1 (See figure on previous page.) Fig. 4 Functionally grouped gene network analysis for Atlantic salmon subjected to Warm & Hypoxic conditions. Network based on a 735 upregulated and b 592 down-regulated differentially expressed genes (DEGs) in salmon exposed to Warm & Hypoxic (WH: 20°C,~70% air saturation) conditions as compared to the Control group (CT: 12°C, 100% air sat.) (n = 6 per treatment, N = 12 total). GO-terms and pathways were obtained through a functional enrichment analysis using the ClueGO plugin in Cytoscape (v3.5.1) [45]. Each node represents a significantly enriched Gene Ontology (GO), KEGG or Reactome pathway (hypergeometric test p < 0.05 with Benjamini-Hochberg correction). The node colour regime visualizes functional groups and processes that share similar genes. The most significant enriched terms of each functional group are illustrated as a summary label. The node circle diameter corresponds to the significance of the enriched pathway (i.e., larger circles correspond to higher significance). The gating and numbering represent the classification of ten functional themes: Heat Shock Response ( # 1), Cellular Stress ( # 2), Oxidative Stress ( # 3), Apoptosis ( # 4), Immune Response ( # 5), Protein Processing & Localization ( # 6), Transcription ( # 7), Proteolysis ( # 8), Catabolic Processes ( # 9) and Cellular Metabolic Processes ( # 10). GO/pathway term annotations were obtained using the GO database for biological process, cellular component, molecular function and immune processes, and the KEGG and Reactome pathway databases. Abbreviations of GO/pathway terms were used for simplification, and the complete list of terms is represented in Additional files 4 and 6    a Linear mixed-effect models (LMEs) were performed for each gene individually to assess the effect of the fixed factor 'treatment' and included the random term  Table 3). Two stress-and hypoxiarelated genes, hif1α (p = 0.042; Fig. 5l) and gstt1 (p = 0.021; Fig. 5m), were only significantly down-regulated in the WH group as compared to the CT group ( Table  3). The gene dnmt1, which is involved in transcriptional regulation, was significantly down-regulated in WN and WH fish in comparison to CT fish (p = 0.007; Fig. 5n, Table 3). Finally, the gene pdk3, which plays a key role in regulating glucose metabolism, was significantly upregulated in the WN group (p = 0.025; Fig. 5o, Table 3), but did not quite reach significance (p = 0.051) in the WH group.

Genes related to apoptosis ( # 4)
Two genes related to apoptotic processes [casp8 (p = 0.023; Fig. 5p) and jund (p = 0.024; Fig. 5q)] were more highly expressed in fish challenged with the WN and WH conditions in comparison the CT group, while jak2 (p = 0.053; Fig. 5r) showed a similar trend ( Table 3). The gene ctsh was significantly down-regulated in WN and WH fish in comparison to CT fish (p = 0.003; Fig. 5s, Table 3).

Genes related to immune response ( # 5)
Three immune-related genes [apod (p = 0.026; Fig.  5t), c1ql2 (p = 0.005; Fig. 5u) and epx (p = 0.019; Fig.  5v)] had higher expression levels in the treatment groups as compared to the CT group. However, these genes had different expression profiles ( Table  3). The gene c1ql2 was significantly up-regulated in both experimental groups as compared to the CT group (Table 3). Whereas the expression of apod was only significantly higher in the WH group, and the expression of epx was only higher in the WN group as compared to the CT group (Table 3). The expression of nckap1l was significantly downregulated in the WN and WH groups as compared to the CT group (p = 0.022; Fig. 5w), and this result was contrary to the up-regulation found for the microarray probe (Additional file 8).
PCA analysis on 24 stress-related ( # 1-4) and 14 immunerelated ( # 5) genes The PCA analysis based on the 24 stress-related genes revealed similar transcript expression profiles for salmon exposed to WN Fig. 6b and Table 1).

Correlation between gene expression and phenotypic traits
The following seven phenotypic characteristics reported in Gamperl et al. [42] were used for correlation analyses with the qPCR gene expression data (RQs): weight, fork length, condition factor (CF), specific growth rate (SGR), hepato-somatic index (HSI), spleen-somatic index (SSI) and relative ventricular mass (RVM).

Correlation analysis for WN fish
In the component map based on the response variables obtained for WN fish (Fig. 6c), the target genes segregated according to their responses into 20 up-regulated and 21 down-regulated genes on opposing directions along Component-1 (CP-1: explaining 37.4% of the variance) (Fig. 6c). The following up-regulated genes contributed significantly to the variance explained by CP-1: 3) (Fig. 6c; Additional file 9). On the contrary, the down-regulated genes that had a significant contribution to the variance of CP-1 were: prdx6 (CP-1 = 11.6), egln2 (CP-1 = 8.  Table 3). The plots are sorted according to six functional themes. Heat Shock Response: a serpinh1, b hsp90aa1, c hsp70; Cellular and Oxidative Stress Response: d hcn1, e cirbp, f calm, g cyp1a1, h egln2, i prdx6, j rraga, k ucp2, l hif1α and m gstt1; Transcription: n dnmt1; Metabolism: o pdk3; Apoptosis: p casp8, q jund, r jak2 and s ctsh; and Immune Response: t apod, u c1ql2, v epx and w nckap1l  Table 1), all of the genes with significant associations to CP-1 are likely to be key drivers of the phenotypic changes in WN fish. Out of the seven phenotypic characteristics for the WN group, only SGR contributed significantly ( Fig. 6c; Additional file 9), and six genes showed a positive correlation (R between 0.5 and 0.65) with SGR (irf2, tapbp, hsp70, hsp90ab1, igfbp2b1 and ndufa1; Fig. 6e).

Discussion
Seasonal fluctuations in water temperature and oxygen content are unavoidable at aquaculture cage-sites and have profound impacts on the growth, survival and welfare of farmed fish [1,3,33,36,40,42]. With global warming, it is predicted that suboptimal temperatures (both high and low) and low water oxygen levels (hypoxia) may become more frequent and severe in coastal areas [7][8][9][10][11], and ultimately more challenging for aquaculture species [46]. However, we still have limited knowledge about the capacity of salmon to tolerate high temperatures in combination with hypoxia, and need to understand how their physiology and immune function are affected. Thus, in this study, Atlantic salmon were subjected to an incremental increase in water temperature (12 → 20°C; at 1°C week − 1 ) under normoxia or moderate hypoxia (~70% of air sat.), conditions that realistically simulate summer conditions in salmon aquaculture sea-cages [33,34,41]. Then, functional genomic approaches were employed (Agilent 44 K salmonid oligonucleotide microarray and Fluidigm Biomark™ qPCR) to assess the hepatic transcriptomic responses to these environmental stressors due to the important role of the liver in numerous biological processes [24,43]. We report that salmon exposed to incremental warming up to 20°C, alone and in combination with moderate hypoxia, showed somewhat similar global gene expression responses that were very distinctive from the control (12°C and normoxic) fish. Overall, we identified a set of 2894 DEPs, out of which 1111 DEPs (38%) were shared between the WH and WN treatment groups. Biological pathway analysis suggested that both treatments increased gene expression with regards to the heat shock response, the unfolded protein response (UPR), endoplasmic reticulum (ER) stress, apoptosis and immune defence. In contrast, a variety of general metabolic processes, proteolysis and oxidation-reduction processes were suppressed. Interestingly, even though the combination of high temperature and moderate hypoxia influenced a number of similar processes, we also identified a unique set of 994 DEPs (34%) that were more strongly dysregulated in WH fish and showed more pronounced impacts on the heat shock response and immune processes. Further, fish exposed to both stressors showed strong correlations between the expression of 19 microarray-identified biomarkers and parameters of growth performance (i.e., weight, length, CF, SGR), which were previously reported to be reduced in these fish [42]. This data reinforces the biological relevance of these genes and pathways, and emphasizes their involvement in phenotypic responses.
Heat shock response, unfolded protein response and endoplasmic reticulum stress In stressful conditions, the expression of highly conserved, ubiquitously distributed, molecular chaperones [Heat Shock Proteins (HSPs)] is initiated to maintain cell function and homeostasis, and to protect tissues and cells from structural damage [47][48][49]. HSPs play a fundamental role in the folding of newly synthesized polypeptide chains, and the refolding and degradation of misfolded proteins to prevent their aggregation and loss of functionality [48,49]. In addition, they are involved in immune system function (e.g., antigen presentation) [50,51], apoptosis [52] and protecting cells from oxidative stress [53]. HSPs are essential regulators of the cellular stress response in aquatic ectotherms [54], and their expression is well characterized in teleost species exposed to thermal stress [13][14][15][16][17][21][22][23] and hypoxia [26,29,30]. After the incremental temperature challenge to 20°C (alone or in combination with hypoxia), we found enriched pathways related to the heat shock response, protein folding and protein stability (Theme # 1) (i.e., 'chaperone-mediated autophagy processes', 'chaperone and protein folding responses', 'STIP1(HOP) binds HSP90 and HSP70:HSP40:nascent protein'). Interestingly, we observed a similar magnitude of up-regulation for genes related to chaperone function in WN and WH fish (i.e., serpinh1, hsp90aa1, hsp90ab1 and hsp70) in the microarray. Of these, serpinh1 (alias hsp47, encoding Serpin H1) was one of the most up-regulated target genes (qPCR validated) in comparison to the control group, and this is in agreement with previous studies on over ten fish species [55]. Serpin H1 binds very specifically to collagens and procollagens to facilitate their assembly and stabilization, and plays an important role in collagen biosynthesis [56]. Moreover, Serpin H1 is involved in the breakdown of reactive oxygen species (ROS) produced during oxidative stress as recently shown in rainbow trout (Oncorhynchus mykiss) [57]. Hence, the increased expression of serpinh1 mRNA in the liver may have assisted in the stabilization of collagen molecules within the extracellular matrix (ECM), and further enabled the elimination of ROS and the maintenance of cellular homeostasis during this thermal challenge. Likewise, the expression of the gene hsp90aa1 (alias hsp90-alpha, encoding Heat Shock Protein 90 alpha) which is the inducible form of HSP90, and its constitutive counterpart hsp90ab1 (alias hsp90-beta, encoding Heat Shock Protein 90 beta), were both upregulated in WN and WH fish as compared to CT fish (although only hsp90aa1 was qPCR validated). These findings are in line with previous studies reporting higher expression of hsp90aa1 mRNA, and minor effects on the constitutive expression of hsp90ab1, after thermal stress in fish [15,58]. HSP90AA1 is an abundant molecular chaperone and is implicated in a wide variety of cellular processes including protection of the proteome, the folding and transport of newly synthesized polypeptides, and the activation of proteolysis of misfolded proteins [54,59,60]. In our study, hsp90aa1 was connected with most of the HSP-related GO/pathway terms (e.g., 'HSF1 activation' and 'dissociation of cytosolic HSF1: HSP90 complex', Figs. 3a and 4a; Additional files 3 and 4), and hence, was an essential component of the unfolded protein response (UPR). Although up-regulation of hsp70 (Heat Shock Protein 70) in the WN and WH groups was not confirmed by qPCR, it was associated with enriched GO/pathway terms such as 'HSP70: HSP40:nascent protein' (Figs. 3a and 4a; Additional files 3 and 4). HSP70 is considered to be a hallmark of the heat shock response, and enhanced transcript expression upon heat stress has been observed in numerous fish species [54,55]. Thus, in the current study, inducible forms of HSP70 may have been involved in molecular chaperone processes within the liver cells of stressexposed fish to assist with the folding of nascent polypeptide chains and the repair and degradation of altered or denatured proteins [61]. In addition, we identified similar enriched pathways in WN and WH fish that are important for ER protein processing (Theme # 6) (i.e., 'translational initiation', 'protein localization/targeting to ER'), and unique UPR pathways in WH fish, that point to the presence of ER-stress. The ER produces a large number of secretory proteins, and has quality-control systems that ensure the correct folding of proteins and their vesicular cellular transport [62]. The ER-stress response initiates the UPR to prevent the assimilation of unfolded proteins and restore ER function [62].
In summary, we found highly active cellular HSP and chaperone-mediated ER-stress responses in the liver of WN and WH exposed salmon that potentially prevented the accumulation of unfolded proteins and maintained cell homeostasis. However, fish exposed to hypoxia in combination with heat stress had larger and more interconnected clusters associated with the HSP and UPR response (Theme # 1) and this indicates that there was a greater induction of these essential cell maintenance processes in the liver of WH fish during this climate change scenario.

Apoptosis
Apoptosis is a process of programmed cell death that allows for the removal of defective cells without the release of intracellular contents, and prevents local tissue inflammation [63]. High temperature and hypoxia can induce oxidative stress that mediates apoptotic processes in fish tissues [20,[64][65][66][67]. Here, we identified similar enriched GO/pathway terms connected to apoptosis (Theme # 4) in WN and WH fish that were associated with up-regulated DEGs (i.e., 'cysteine-type endopeptidase activity involved in apoptotic signaling pathway' and 'regulation of nitrosative stress-induced intrinsic apoptotic signaling pathway'). However, the WN fish had more enriched GO/pathways terms related to apoptotic processes, and these had a higher interconnectivity in the functional network (~92 terms) as compared to WH fish (~10 terms). In addition, WN fish had enriched heat shock response (e.g., unfolded protein binding) and immune (e.g., MHC protein complex binding) pathways that were linked with apoptotic processes. This suggests that enhanced HSP transcript expression (i.e., hspaa1 and hspd8) may have been associated with cell resistance to apoptotic cell death [52].
Nevertheless, in both warming groups, we found several differentially expressed genes from the microarray experiment that were associated with these enriched apoptosis-related pathways (i.e., casp8, jund and jak2). Amongst them, the gene casp8 was the most frequently occurring transcript in enriched apoptosis-related GO/ pathways (Additional files 3 and 4), and was significantly up-regulated in WN and WH fish as compared to CT fish (and qPCR validated). CASP8 initiates a protease cascade that induces receptor-mediated extrinsic cell death [68]. On the contrary, the significant upregulation of the transcript jund (encoding Transcription factor JunD) in WN and WH fish indicates that antiapoptotic processes were also activated to protect cells from senescence or apoptosis by acting as a modulator of the p53 signalling cascade [69].
Collectively, these findings suggest that several processes were involved in initiating the apoptosis of cells potentially damaged by oxidative stress in the liver of WN and WH fish, but also that anti-apoptotic pathways likely prevented extensive programmed cell death which could have resulted in hepatic necrosis or impacted liver function [1,63,66,67]. Interestingly, moderate hypoxia combined with high temperature resulted in a less pronounced activation of apoptotic processes (Theme # 4) as compared to fish that experienced high temperature alone. At present, we do not have an explanation for this finding.

Immune response
Several studies have shown that exposure to elevated temperatures [15,[17][18][19]22] or hypoxia [6,70] affects immune-related gene expression in fish. However, these studies focused on the effects of acute or chronic temperature increases, or hypoxia, as individual stressors. Here we report that an incremental increase in temperature under normoxia or moderate hypoxia resulted in the up-regulation of several genes linked to the innate and adaptive immune system (Theme # 5) such as apod, c1ql2, casp8, epx, mhcii, il8, tapbp and tnfrsf6b in the microarray (with apod, c1ql2 and casp8 qPCR validated). Interestingly, salmon exposed to moderate hypoxia and high temperature had more up-regulated genes that were associated with enriched GO/pathway terms of the innate immune response (i.e., 'neutrophil and granulocyte chemotaxis and granulation', 'IL-17, IL-4 and IL-13 signalling'), and antiviral responses (i.e., 'herpes simplex virus 1 infection'). Further, WH fish displayed a more distinct immune-related gene expression profile as compared to CT and WN fish. This shift could be attributed to a stronger differential expression of several genes in WH fish such as apod, camp-a, c1ql2, il8, mhcii and tnfrsf6b, and some of these correlated strongly with health-related parameters. For instance, apod (encoding Apolipoprotein D) was more up-regulated in WH fish as compared to WN fish (WN: 4.4-fold, WH: 7.7-fold; qPCR validated), and was positively correlated with spleen-somatic index. The extracellular glycoprotein APOD has multiple functions that involve the immune response, chemoreception, proteolysis and lipid oxidation [71], and potentially plays a fundamental role in cellular-stress and immune responses. The higher expression of c1ql2 (encoding C1q-like Protein 2) transcripts (WN: 6.4-fold, WH: 12.3-fold; qPCR validated) suggests that the classical complement system pathway was also activated in these treatment groups. This is an essential innate defence mechanism in teleost fish that detects and destroys invading pathogens by bacteriolysis [72]. For example, acclimation to 20°C for 57 days increased the lytic activity of the total complement system in the plasma of rainbow trout [73]. Interestingly, the gene c3 (encoding Complement C3), an important component of the alternative complement pathway, was not significantly affected by exposure to elevated temperatures in the qPCR assessment. Consistent with this finding we found that there was no change in plasma hemolytic activity of the alternative pathway in the same WN and WH fish 24 h after an intraperitoneal injection with bacterial and viral antigens in comparison to CT fish [74]. Oku and colleagues [19] reported that the classical pathway of the complement system was partially activated in the liver of rainbow trout reared at 22°C vs 14°C whereas there was no activation of the alternative pathway or induction of final cytotoxic activity. The classical pathway links innate and adaptive immunity as C1q has the ability to bind to aggregated IgG or IgM on immune complexes [75]. Indeed, we found enrichment of 'MHC complex binding' processes and a trend for higher expression of mhcii (alias hla-a, encoding Major Histocompatibility Complex II) in WH fish (WN, 1.0; WH, 1.5-fold; but not qPCR validated). Thus, the applied thermal challenge may have activated adaptive immune responses, on which fish rely more at high temperatures [76]. The GO/pathway term network also indicated that several HSP transcripts were associated with 'antigen processing and presentation' (e.g., hsp90aa1, hsp90ab1, hspa8) and with 'MHC class II binding complex' (e.g., hsp90aa1 and hspa8). HSPs have been shown to act as chaperones for cytosolic peptides involved in MHCdriven antigen presentation to T-lymphocytes [50]. For instance, the cytosolic chaperone HSP90 is associated with peptide binding to MHC-class I and MHC-class II [77,78]. Consequently, given the observed up-regulation of hsp90aa1 transcripts in WN and WH fish as compared to CT fish, and this HSP's relation to 'antigen processing pathways', it is possible that MHC peptide complex assembly for antibody production was enhanced [50,51,77,78].
Finally, WN and WH fish had up-regulated genes connected to collagen/ECM degradation (i.e., 'activation of matrix metalloproteinases' or 'collagen degradation by MMPs') and a trend for higher expression of mmp9 (encoding Matrix Metalloproteinase 9) (WN, 2.0; WH: 1.5-fold; not qPCR validated). This suggests that these environmental conditions stimulated tissue remodeling processes [79]. Matrix metalloproteinases (MMPs) are endopeptidases that cleave all structural elements of the ECM and are responsible for physiological and pathophysiological tissue remodelling [79,80]. As such pathological remodeling processes in damaged liver tissues due to apoptosis were likely activated.
In summary, our results suggest that the constitutive expression of immune-related genes was induced to either prepare for more numerous and/or virulent pathogens at warm temperatures [81,82] as a 'pre-adaptation', or to initiate immune defence responses against invading pathogens or existing infections. Further, they show that the combined stressors of high temperature and moderate hypoxia had a greater impact on hepatic immunerelevant transcript expression. However, no significant clinical signs of infection or mortalities were recorded in this experiment [42], and when salmon from both warming scenarios were held at 20°C for 4 weeks (WN and WH groups) and challenged with a multivalent vaccine (Forte V II; containing both bacterial and viral antigens), their capacity to mount an innate immune response was not impaired and they reached a similar magnitude of antibacterial immune-related gene expression as compared to CT fish [74]. Nonetheless, it is clear that increasing temperatures due to climate change [9] will become more challenging for Atlantic salmon. For example, plasma cortisol levels, which are known to modulate and suppress immune function at higher concentrations [83,84], were significantly increased (tõ 30 to 40 ng mL − 1 ) in Atlantic salmon exposed to an incremental temperature elevation to 22°C (Zanuzzo FS, Peroni EFC, Sandrelli RM, Beemelmanns A, Dixon B, Gamperl AK. The impacts of increasing temperature and moderate hypoxia on the stress physiology of Atlantic salmon (Salmo salar) (In prep.)), with approx. 15% mortality at this temperature and 33% when temperatures reached 23°C [42]. Hence, temperature elevations above 20°C will likely have a stronger impact on physiological stress and immune competence, and ultimately the disease resistance of Atlantic salmon. Clearly, additional research using live pathogen exposures is needed to determine whether the susceptibility of salmon to infectious diseases is impacted by warmer and more hypoxic environments.

Oxidative stress
When ectotherms are exposed to warming (hyperthermia) their mitochondrial respiration is increased, and this can result in accelerated mitochondrial ROS formation, oxidative stress and cellular damage [64,85,86]. Cellular oxidative stress due to prolonged hyperthermia can cause impaired mitochondrial bioenergetics and structural alterations of cells and tissues [67,87].
Tolerance to cellular oxidative stress is provided by an effective antioxidant system [86,88], as well as an HSP response which attempts to maintain protein folding and mitochondrial integrity, and support cell function and survival [52,53]. In our study, WN fish appeared to have an enhanced oxidative stress response (Theme # 3) (i.e., 'regulation of response to stress' and 'response to oxidative stress') while WH fish showed up-regulated redoxrelated pathways (i.e., 'antioxidant activity', 'positive regulation of release of cytochrome c from mitochondria' and 'thioredoxin reduction'). Thus, the activation of these particular oxidative stress responses, in addition to activated HSP, UPR, and ER-stress responses, appear to be critical in preventing cellular damage and maintaining cell homeostasis during warming.
Surprisingly, all of the measured oxidative stress and hypoxia-sensitive target genes that were selected for qPCR validation (i.e., cirbp, calm, cyp1a1, egln2, prdx6, rraga, ucp2) were down-regulated in the liver of WN and WH fish at 20°C as compared to CT fish. Further, some GO/pathway terms connected to cellular oxidative stress (Theme # 3) (i.e., 'oxidoreductase activity') were associated with down-regulated genes in both warmed groups. Olsvik et al. [14] also reported that the expression of several genes encoding for proteins with an oxidative stress-protective and/or hypoxia sensing function (i.e., sod1, gr, cyp1a1, hif1α) was significantly reduced in the liver of Atlantic salmon exposed to prolonged high temperature stress (19°C vs 13°C for 45 days). In accordance, we also found that the gene hif1α (alias hif-1a, encoding Hypoxia-Inducible Factor 1 alpha) was significantly down-regulated in WH fish (by 0.57-fold, qPCR validated), while it fell just short of being significant for WN fish (WN: 0.60-fold; p = 0.054). The gene hif1α encodes for a master hypoxia-responsive transcriptional regulator involved in various cellular processes such as energy metabolism, apoptosis, proliferation and increased oxygen delivery [24,89,90]. The expression of hif1α is considered to be an evolutionarily conserved hypoxic biomarker due to its up-regulation after acute hypoxia (i.e., hours) in several fish species such as Eurasian perch (Perca fluviatilis) [91], Atlantic croaker (Micropogonias undulatus) [92] and zebrafish (Danio rerio) [28]. However, the data are not as consistent regarding the effects of prolonged hypoxia. For example, while chronic hypoxia induced the up-regulation of hif1α transcripts in the ovary of Atlantic croaker (21 days at 55% DO) [92] and in the liver of sea bass (Dicentrarchus labrax) (15 days at 51% DO) [93], the expression of this transcript was down-regulated in the muscle and not affected in the liver of perch (15 days at 30% DO) [91]. Further, while this gene's expression was not altered in the liver of Atlantic salmon exposed to 4-5 mg O 2 L − 1 for 120 days at 12°C, long-term exposure to 17°C and 19°C (45 days) as compared to 13°C resulted in a lower expression of hif1α in the liver of Atlantic salmon [14]. Moreover, Heise [86] showed that while the DNA binding activity of the transcription factor HIF-1 in North Sea eelpout (Z. viviparous) liver cells was elevated during mild heat exposure (18°C), its function appeared to be impaired when this species was exposed to more severe temperature stress (22-26°C). This author hypothesized that a more oxidized redox state during very high temperatures could interfere (i.e., 'switch-off') with the hypoxic signalling response, and thus, prevent the complex HIF-1 induced physiological response. This hypothesis may be supported by the significant down-regulation of egln2 (alias phd1, encoding Egl-9 Family Hypoxia Inducible Factor 2) observed in both warmed groups (qPCR validated) in this study, as it encodes for cellular oxygen sensing enzymes responsible for the post-translational regulation of HIF-1α proteins [94][95][96]. Three EGL-Nine homologs (EGLN1-3) regulate the abundance of HIF-1α proteins through proline hydroxylation and consequent proteasomal degradation [96], and were shown to be involved in signalling responses in the brain of large yellow croaker (Larimichthys crocea) exposed to acute hypoxia [97]. The effects of temperature stress on egl-9 homolog transcript expression have not previously been reported. Interestingly, we found a down-regulation of egln2 transcripts in both the WN and WH groups, and this may imply that there was temperature-dependent posttranslational regulation of HIF1α in the liver of our Atlantic salmon. In addition, the down-regulation of calm (alias cam, encoding Calmodulin) in WN and WH fish as compared to CT fish (qPCR validated) suggests that there might have been Ca 2+ /calmodulin kinasedependent transcriptional regulation of HIF-1 [98]. The Ca 2+ /calmodulin pathway was supressed in the liver of hypoxia tolerant gynogenetic blunt snout bream (Megalobrama amblycephala) [31], and thus, may have been important in mediating hypoxia acclimation in salmon. Based on these findings, it is possible that the lower expression of hypoxia sensitive genes (hif1α, calm, egln2) in the liver of WH fish may have been caused by: i) the moderate level of hypoxia (~70% air saturation); ii) an acclimation response to prolonged~8-10 weeks of hypoxia/temperature stress; and/or iii) a negative feedback loop due to the accumulation of Hif1α proteins. However, further research is needed to gain a better understanding about how the HIF1α pathway in Atlantic salmon is modulated by these two important environmental stressors.
In relation to the above findings, the gene cyp1a1 which encodes for Cytochrome P450 1A1 was downregulated in fish from the WN and WH treatments as compared to CT fish (qPCR validated). Hypoxia and temperature stress can alter transcript levels of CYP1A, which is involved in the oxidation of many substrates and considered to be a vital molecular biomarker for various stressors in the aquatic environment [99]. In line with our results, elevated temperatures (17-18°C) caused a down-regulation of cyp1a mRNA in the liver of Atlantic salmon [14], and chronic hypoxia decreased the expression of this gene in Atlantic cod (six weeks at 46% O 2 saturation) [100] and Atlantic croaker (four weeks at 1.7 mg DO L − 1 ) [101]. The expression of ucp2 was also downregulated in WN and WH fish (qPCR validated) and this gene codes for Mitochondrial Uncoupling Protein 2 (UCP2) which uncouples oxidative phosphorylation from ATP synthesis resulting in energy dissipation [102]. Decreased transcript levels of ucp2 with increasing temperatures (15-25°C) were previously reported in the gill and liver of pikeperch (Sander lucioperca) [58], and suggest that UCP2 has a thermogenic function [102]. Further, gilthead sea bream (Sparus aurata) showed a significant decrease in the expression of ucp2 in the whole blood after acute (1 h) exposure to 18-20% oxygen saturation [103]. Under stressful conditions, reduced UCP mediated uncoupling (respiration uncoupling) may result in the attenuation of mitochondrial ROS production, and a lower expression of ucp2 could have been part of a feedback-induced decrease in ROS synthesis for cell protection [104]. Indeed, Gerber and co-workers reported that Atlantic salmon acclimated to 20°C had reduced cardiac mitochondrial ROS production in comparison to fish acclimated to 12°C [105]. Thus, alterations in mitochondrial function at high temperatures may be an important mechanism for thermal acclimation and thermal tolerance in this species.
Another highly down-regulated gene in WN and WH fish was prdx6 (qPCR validated), which encodes for Peroxiredoxin-6. This protein is important for phospholipid homeostasis, lipid peroxidation repair, and inflammatory signalling [106]. Up-regulation of prdx6 upon heat stress to protect the cell from oxidative stress has been reported in other marine animals [107,108]. However, while Antarctic emerald rockcod (Trematomus bernacchii) exposed to warming temperatures had slightly increased expression of the prdx6-b paralog in the liver, expression of the prdx6-a paralog was down-regulated [108]. This suggests that prdx6 paralogs may respond differently to temperature changes, and that the transcriptional responses of this gene and its paralogs upon warming and exposure to hypoxia deserve further investigation.
Finally, the expression of cirbp (which encodes for Cold-Inducible RNA-Binding Protein) was downregulated to a similar extent in WN and WH fish (qPCR validated), and this transcript was connected with many GO/pathway terms including 'mRNA stability' and 'mRNA catabolic process'. The expression of cirbp has been reported to be up-regulated upon cold water exposure, and mild hypoxia, while it is decreased in response to heat stress and chronic hypoxia in vertebrates [55,109]. The cold-shock protein CIRBP acts as mRNA chaperone, is implicated in multiple cellular processes (i.e., cell proliferation, survival and apoptosis), and is considered to be a general stress-response protein affected by temperature, hypoxia and UV radiation [109]. The lower expression of cirbp in this study is in accordance with several heatstress studies on salmonid fishes [55].
Taken together, the up-regulation and enrichment of pathways related to oxidative stress (Theme # 3), HSPresponse (Theme # 1) and apoptosis (Theme # 4) in WN and WH fish suggest that the induction of antioxidant enzymes and redox pathways was an important defense mechanism in these fish. However, the simultaneous down-regulation of several key genes related to the oxidative stress response (Theme # 3) also suggests a potential decrease in its effectiveness when salmon are exposed to a slow incremental temperature increase to 20°C.

Cellular metabolism
The abiotic factors of temperature and water oxygen level have a profound influence on the allocation of energy to maintenance versus growth in fish [3,38,110]. A reduction in metabolic processes can conserve energy during stressful conditions as imposed by thermal challenges [14,17,20,22,32,65,111,112] or hypoxia [25,28,113,114]. In this study, the expression of genes related to a variety of highly interconnected cellular metabolic processes was suppressed in WN and WH salmon. Amongst many others, the GO/pathway terms of down-regulated genes were connected to aerobic respiration [i.e., 'tricarboxylic acid cycle (TCA)'], carbohydrate metabolic process (i.e., 'glucose 6-phosphate metabolic process'), and small-molecule metabolic process (i.e., 'organic substance biosynthetic process', 'fatty acid catabolic process', 'lipid metabolic process'). The downregulation of genes associated with the TCA cycle in the liver of WN and WH fish may indicate a shift from aerobic oxidation to anaerobic glycolysis, as was shown in the Tambaqui (Colossoma macroppomum) exposed to predicted IPCC climate scenarios [32]. The decrease in the expression of hepatic gck (detected in the microarray), which encodes for the enzyme Glucokinase, suggests that there was also a reduction in glycolytic processes in the liver of WN and WH fish. This result is in agreement with the response of Tambaqui exposed to extreme climate scenarios [112]. In contrast, expression of the gene pdk3 was up-regulated in WN fish (qPCR validated). This gene codes for Pyruvate Dehydrogenase Kinase 3 (PDK3), which acts together with PDK1, PDK2 and PDK4 isoenzymes to regulate glycolysis and glucose homeostasis during starvation [115].
In our study, WH salmon had reduced food consumption, and a lower feed conversion ratio and growth, as compared to fish of the WN and CT groups [42]. The reduced feed intake and feed conversion efficiency at high temperatures could have had a major impact on the redistribution of energy stores and the amount of glycogen and lipids stored in the liver. Temperature modulates lipid metabolism, and stored lipids in the liver are increasingly used for the maintenance of energy metabolism during thermal stress [38,110]. For example, Atlantic salmon reared at 17-19°C for 45 days showed reduced liver lipid and triacylglycerol (TAG) stores as compared to fish maintained at 13°C, and this suggests that the reallocation and/or depletion of endogenous lipid stores occurs during prolonged high temperature exposure [38]. Furthermore, Atlantic salmon held at 18°C vs 12°C for 1 month showed a decline in plasma amino acids (glutamine, tyrosine and phenylalanine) and a decreased lipid status (unsaturated fatty acids, lipids and phospholipids), suggesting that energy stores were mobilized [110]. In the current study, the expression of genes associated with 'fatty acid and lipid metabolic processes' was lower in WH fish as compared to CT fish, and thus, it appears that long-term exposure to high temperatures and hypoxia may result in reduced lipid and fatty acid biosynthesis. This hypothesis would be consistent with a recent study on rainbow trout exposed to an incremental temperature increase to 24°C. These fish showed a temperature-dependent down-regulation of hepatic genes related to energy metabolism [20]. Thus, a down-regulation of pathways related to carbohydrate, protein and fatty acid metabolism in the liver of stressed salmon may reflect the suppression of metabolic processes, and agrees with the important role played by the liver in cellular metabolism and biosynthetic activities in fishes [24,38,43].
During hypoxia, metabolic responses to ensure cell survival involve readjustments that decrease ATP demands to match the reduced capacity for ATP production [24]. Moreover, prolonged temperature stress and low oxygen reduce protein synthesis, and this leads to reduced growth and metabolic depression in Atlantic salmon [14]. These findings are consistent with our results. We found positive correlations between the decreased expression of 12 down-regulated genes (e.g., cirbp, calm, egln2, hif1α, ucp2, gstt1, prdx6, rraga, etc.) and reduced growth performance (i.e., weight, length and SGR) predominantly in WH fish [42], and this association highlights the biological relevance of these hepatic transcriptional responses.
Collectively, these findings suggest that the combination of high temperature stress and moderate hypoxia resulted in transcriptional responses in the liver that contributed to metabolic suppression in our salmon. This metabolic suppression may have been at least partially needed to balance the energetically costly processes that were invoked to maintain cell homeostasis (i.e., HSP, UPR, ER-stress and apoptosis).

Transcriptional regulation and epigenetic mechanisms
Temperature stress in the WH and WN groups induced a similar down-regulation of the gene dnmt1 as compared to CT fish (qPCR validated), and this gene codes for DNA (cytosine-5)-methyltransferase, an enzyme essential for maintaining DNA methylation marks after mitosis [116]. DNA methylation is an important epigenetic regulatory mechanism of transcription (Theme # 7), and down-regulation of dnmt1 indirectly suggests that genome-wide changes in DNA methylation levels may have been involved in regulating these large-scale gene expression responses (i.e.,~2894 DEPs). Indeed, Beemelmanns et al. [117] found that the same treatments (i.e., WH and WH) affected the methylation of CpG sites of the microarray-identified genes related to temperature stress (serpinh1, cribp), oxidative stress (prdx6, ucp2), apoptosis (jund) and metabolism (pdk3). Several of these changes in CpG methylation were highly correlated with the transcript expression changes reported here, and thus, reinforce their importance as 'epimarkers' that regulate transcription upon temperature and hypoxic stress in Atlantic salmon [117].

Conclusions
We identified numerous transcriptional changes (i.e., 2894 DEPs) in the liver of Atlantic salmon exposed to an incremental temperature increase (12 → 20°C; at 1°C week − 1 ) alone or combined with moderate hypoxia (7 0% of air sat.); the latter simulating summer conditions in salmon aquaculture sea-cages [33,34,41]. Both these treatments induced biological processes related to the maintenance of cellular homeostasis (i.e., HSP, UPR, ERstress response) and the apoptosis of damaged cells, and stimulated immune gene expression (both innate and adaptive); but also compromised the oxidative stress response and led to a reduction in the expression of a variety of genes related to metabolic and proteolytic processes. Importantly, we identified a unique set of 994 DEPs (34%) that were strongly dysregulated in WH fish, and showed that this condition had a more distinct and pronounced impact on the heat shock response and immune processes (e.g., more up-regulated genes that were associated with enriched GO/pathway terms of the innate and antiviral immune responses) as compared to that seen for WN fish.
Transcriptomic techniques allow for the highthroughput identification of genes that are sensitive to particular conditions, and can be used as biomarkers for the detection and quantification of stress levels and stress tolerance. In this context, the development of diagnostic biomarkers for quantifying the impact of environmental stressors on an organism's physiology and health has received increased attention by the aquaculture industry and for use in ecological surveys [55,58,118]. Our results agree with recent studies which show that the genes serpinh1, hsp90aa1 and cirbp are reliable molecular biomarkers for the detection and quantification of thermal stress in salmonids [55,58,119]. However, we found that the expression of hsp70 was very variable between individuals, and in accordance with previous findings, should not be considered as a stress biomarker alone [120].
We also report that the marked differential expression of 19 microarray-identified genes upon high temperature and hypoxia exposure showed strong associations with important phenotypic characteristics. Hence, these genes may not only be useful as molecular biomarkers of thermal stress, but also as candidate genes for the development of thermal phenotype-relevant genetic markers [e.g., single nucleotide polymorphisms (SNPs) for marker-assisted selection of heat stress resistant broodstock], protein-based diagnostic assays [e.g., Enzymelinked Immunosorbent Assay (ELISA)] and for the detection of epigenetic markers ('epimarkers') that can predict thermotolerance [117].

Methods
This experiment was performed as part of the 'Mitigating the Impacts of Climate-Related Challenges on Salmon Aquaculture (MICCSA)' project, and a detailed description of the experimental protocol and of the data on growth characteristics and mortality are published in Gamperl et al. [42]. All experimental procedures described herein were approved by the Institutional Animal Care Committee of Memorial University (Protocol # 16-90-KG) and followed guidelines of the Canadian Council on Animal Care. All sections of this study adhere to the ARRIVE Guidelines for reporting animal research [121].

Animal husbandry
The experiment was performed from March to August 2017 at the Laboratory for Atlantic Salmon and Climate Change Research (LASCCR), Memorial University, St. John's, NL, Canada. Post-smolt Atlantic salmon (~1.5 years old) of Saint John River (NB, Canada) origin obtained from Northern Harvest Sea Farms Ltd. were implanted with Passive Integrated Transponder (PIT) tags (Loligo® Systems ISO 11784 certified, Viborg, Denmark), then randomly distributed into six 2.2 m 3 circular indoor fiberglass tanks receiving seawater (32 ppt salinity) at 15 L min − 1 . The fish were initially acclimated for four weeks under optimal conditions (~100-110% air saturation, 12°C, 32 ppt salinity, 14 h light: 10 h dark photoperiod) and fed a ration of 1% body weight day − 1 with a commercial salmon feed (5 mm, EWOS Dynamic S, EWOS Canada Ltd., Surrey, BC, Canada).

Experimental protocol
Two tanks with 60 fish per tank (average mass 137.6 ± 1.3 g; mean ± S.E.), were randomly assigned to each of three groups as shown in Fig. 1: (1) CT, constant temperature of 12°C and~100-110% air sat. for the duration of the experiment; (2) WN, incremental temperature increase (12 → 20°C at 1°C week − 1 ) at1 00-110% air sat.; and (3) WH, decrease in oxygen content to~70% air sat. (daily range~65-75%) over one week, followed by two weeks of acclimation to this oxygen level, and then the same temperature regimen at7 0% air sat. The weekly temperature increases in the tanks of the WN and WH treatment groups were 0.3°C (from days 1 to 3), 0.1°C on day 4, and then no change from days 5-7. The temperature and dissolved oxygen level in the tanks were monitored daily (YSI, ProODO, Yellow Springs, OH, USA), and ammonia and nitrite levels in the tanks were measured weekly (LaMotte test kits, Chestertown, MD, USA). During the experiment, the salmon were carefully fed by hand to satiation twice daily (at 9:00 and 15:00) with the same commercial salmon feed. See Gamperl et al. [42] for more specific information on these experimental protocols.
In the current study, we sampled eight salmon per treatment group after the fish were exposed to 20°C for three days (four fish per tank replicate, n = 8 per treatment, N = 24 fish total). The number of 24 samples was considered sufficient to achieve statistical robustness and power (80%) to detect a significant effect (p < 0.05) with an estimated medium-large effect size (f 2 = 0.43) [122]. All the fish used in this study were euthanized in aerated seawater containing 0.4 g L − 1 of MS-222 (tricaine methanesulphonate; Syndel Laboratories, Nanaimo, BC, Canada) followed by cerebral concussion. Fish were dissected and 200-300 mg of liver tissue was collected, flash-frozen in liquid nitrogen and then stored at − 80°C. The fish's weight [g], fork length [cm], liver mass [g], spleen mass [g] and ventricle mass [g] were measured, and CF, SGR, HSI, SSI and RVM indices were calculated. These growth and physiological measurements are reported in Gamperl et al. [42] and were used in the current study to perform correlation analyses with the gene expression data.

RNA extraction, DNase treatment and column purification
For RNA extraction, 100 mg of liver tissue was disrupted and homogenized in 800 μL of QIAzol-Lysis Reagent (QIAGEN, Germantown, MD, USA) for 2 min at 20 Hz using a TissueLyzerII with 5 mm stainless steel beads (QIAGEN, Mississauga, ON, Canada) according to the manufacturer's instructions. To remove lipid/protein contamination, we performed an additional precipitation step according to the protocol of Xu et al. [123] with minor changes. Briefly, crude RNA samples (100 μg) were mixed with an equal volume of Acid-Phenol: Chloroform (5:1 solution, pH 4.5, Ambion/Life Technologies, Waltham, MA, USA) and centrifuged at 17, 000 g and 7°C for 20 min. Then, 190 μL of the aqueous phase was precipitated with 0.1 volumes of 3 M sodium acetate (pH 5.5; Invitrogen/Thermo Fisher, Vilnius, Lithuania) and 2.2 volumes of ice-cold 100% ethanol (Greenfield Global, Brampton, ON, Canada) at − 80°C for 12 h, and centrifuged at 17,000 g and 7°C for 30 min. The RNA pellets were then washed in 1 ml of 70% ethanol, centrifuged at 17,000 g and 7°C for 20 min, airdried at room temperature for 10 min, and dissolved in 100 μL of nuclease-free water (Invitrogen/Life Technologies, Grand Island, NY, USA) at 55°C for 5 min. To remove any remaining genomic DNA, 40 μg of precipitated RNA was incubated for 10 min with 6.8 Kunitz units of DNase I (RNase-Free DNase Set, QIAG EN, Mississauga, ON, Canada) and 1 × of the manufacturer's buffer. DNase-treated RNA samples were then column-purified using the RNeasy Mini Kit (QIAGEN, Hilden, Germany) according to the manufacturer's guidelines. RNA integrity was tested with gel electrophoresis (1% agarose) and RNA purity and yields were measured by A260/280 and A260/230 NanoDrop UV spectrophotometry (NanoDrop, Wilmington, DE, USA). All column-purified samples had A260/280 ratios between 2.0 and 2.2 and A260/230 ratios between 1.9 and 2.3.

Microarray protocol Microarray experimental design and hybridization
Six individual samples from each treatment (three per replicate tank; n = 6, N = 18 total) were selected to perform a transcriptome microarray study using a common reference design. This sample size number was based on prior microarray studies by our group [70,[124][125][126], and sample-size and power calculations for microarrays using the sizepower function of the Bioconductor package in the R-software [127]. Given the experimental design, we estimated that the power of our statistical analysis was 88%, and that we had the ability to detect2 000 DEGs. An Agilent 44 K Atlantic salmon oligonucleotide array platform (GPL11299; Agilent Technologies, Mississauga, Canada), developed by the consortium for Genomic Research on All Salmonids Project (cGRASP) [44], was used as described in Xue et al. [126] according to the MIAME guidelines [128]. Anti-sense amplified RNA (aRNA) was in vitro transcribed from 1 μg of column-purified RNA using Ambion's Amino Allyl MessageAmp II aRNA Amplification kit (Invitrogen/Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. Thereafter, a common reference pool was generated by combining aRNA from all 18 samples (10 μg from each sample), and 20 μg of aRNA per individual sample or reference was precipitated and re-suspended in coupling buffer. The resulting aRNA was labelled with either Cy5 (individual samples) or with Cy3 (common reference) (GE Health-Care, Mississauga, ON, Canada) through a dye-coupling reaction, following the manufacturer's instructions. The labelling efficiency and concentration of aRNA were determined using spectrophotometry (microarray function of NanoDrop), and the efficiencies ranged between 40 and 60 dye molecules per 1000 nt for all samples. Thereafter, 825 ng of labelled aRNA per sample was mixed with an equal amount of labelled common reference aRNA, and the resulting pools were fragmented and co-hybridized to 44 K microarrays at 65°C for 17 h with 10 rpm rotation using an Agilent hybridization oven (Agilent, Mississauga, ON, Canada). After incubation, the array slides were washed according to the manufacturer's instructions and dried through centrifugation at 200 g for 5 min at room temperature.

Microarray data acquisition
Each 4 × 44 K Agilent microarray slide was scanned at 5 μm resolution with 90% of laser power using a ScanArray Gx Plus scanner and ScanArray Express software (v4.0; Perkin Elmer, Woodbridge, ON, Canada). The Cy3 and Cy5 channel photomultiplier tube (PMT) settings were manually adjusted to balance the fluorescence signal between channels and the four arrays on each slide. The extraction of the resulting raw fluorescence intensity data stored in the TIFF images was performed with Imagene 9.0 (BioDiscovery, El Segundo, CA, USA). Quality controls (background correction and removal of low-quality/ flagged spots) and data transformation (log 2 -transformation and Loess normalization) were performed in R using mArray from the Bioconductor package [129], and with scripts adapted from Booman et al. [130] and Xue et al. [126]. Microarray probes that were absent (non-detected) in more than 25% of the arrays were omitted, and missing values from these undetected probes were imputed using the LSimpute package in R [131].

Microarray data analysis
To identify significantly up-or down-regulated genes in response to the WN or WH treatments as compared to the CT, we applied a permutational SAM test with a false discovery rate (FDR) of < 5% [132] using the implemented siggenes function of the Bioconductor package in R [133]. The SAM-identified DEPs were re-annotated using the contiguous sequences (contigs) which contain the 60mer oligonucleotide probe sequences [44]. We performed BLASTx searches of contig sequences against the Swiss-Prot database and BLASTn searches of 60mer probes against the NCBI nr/nt database with E-value < 1 × 1e-5 as filter criteria as described in Umasuthan et al. [124]. The probe annotations were also updated by homology sequence searches in genome annotation databases for S. salar and O. mykiss, and the assignment of gene symbols for the microarray probes was performed with HUGO Gene Nomenclature Committee (HGNC; https://www.genenames.org/) and/or GeneCards (https://www.genecards.org/) databases according to Umasuthan et al. [124].
A hierarchically clustered heat map was generated with normalized log 2 ratios for 2894 SAM-identified DEPs (FDR < 5%) using the pheatmap R package [134] and applying Pearson correlations and Ward's agglomerative linkage method (ward.D2) as a cluster algorithm.
Functionally organized GO/pathway term network analysis GO/pathway term enrichment and network analyses were conducted for up-and down-regulated DEGs from the WN or WH probe lists (i.e., gene symbol annotation of SAM identified DEPs) using the ClueGO plugin in Cytoscape (v3.5.1) [45]. However, only non-redundant and annotated DEGs were recognized and included (WN group: 377 up-regulated DEGs and 798 downregulated DEGs; WH group: 735 up-regulated DEGs and 592 down-regulated DEGs). To identify enriched gene clusters related to the WN and WH treatments, an enrichment/depletion analysis was performed using a twosided hypergeometric test and the Benjamini-Hochberg p-value correction method [135]. The genes were mapped to GO-term databases for biological, cellular, molecular and immune processes [136] as well as KEGG [137] and Reactome pathway databases [138] in August 2019. Four functionally organized GO/pathway term networks were created with Kappa-statistics (threshold of 0.4), GO-term fusion strategy, and a medium specificity level to minimize the complexity. Due to discrepancies in complexity between the networks, we selected a p-value cut-off of < 0.05 for less complex networks (i.e., up-regulated DEGs), while a lower cut-off of p < 0.005 was applied for highly complex networks (i.e., downregulated DEGs). The enriched terms for each functional group were ranked based on their significance level, and the most significant terms were illustrated as a summary label in the produced networks. The determined GO/ pathway term groups were categorized into the following functional themes: Heat Shock Response ( # 1), Cellular Stress ( # 2), Oxidative Stress ( # 3), Apoptosis ( # 4), Immune Response ( # 5), Protein Processing & Localization ( # 6), Transcription ( # 7), Proteolysis ( # 8), Catabolic Processes ( # 9) and Cellular Metabolic Processes ( # 10). Finally, the non-redundant significantly enriched GO-terms of biological processes were summarised using the REVIGO cluster algorithm [139] and visualized with dot plots.

qPCR protocol Gene selection and primer design
The selection of 41 microarray-identified GOIs for validation purposes included significantly up-and downregulated DEGs that were associated with important functional themes determined by the aforementioned GO/ pathway term network analysis ( # 1-5, # 7, # 10; Table 2). We designed paralog-specific qPCR primers according to specifications detailed in Caballero-Solares et al. [125] using the Primer3web platform (v4.1.0; http://bioinfo.ut. ee/primer3/). Details on qPCR primer sequences, accession numbers, amplicon sizes and amplification efficiencies are represented in Additional file 10. To identify existing paralogs, and to verify the identity of each GOI, we performed BLASTn searches against the nonredundant nucleotide (nr/nt) and the expressed sequence tag (EST) databases of NCBI (year: 2018) (Additional file 11). For nine GOIs (camp-a, cat, cyp1a1, epx, gck, hif1α, hsp70, il8 and mmp9) the paralog-specific primer sequences were obtained from previous studies [140,141] or the Genome Canada Funded Genomic Applications Partnership Program (GAPP) project # 6604 (Biomarker Platform For Commercial Aquaculture Feed Development) database, given that they were targeting the identical probe sequences (Additional files 10 and 11).

cDNA synthesis and primer efficiencies
First-strand cDNA templates for qPCR were synthesized in 20 μL reactions from 1 μg of DNaseI-treated, columnpurified, total RNA with the QuantiTect®Reverse -Transcription Kit (QIAGEN, Mississauga, ON, Canada) following the manufacturer's protocol. qPCR primer quality testing was conducted according to previously published protocols [125,126,130,142]. For the evaluation of primer performance and amplification efficiencies for each primer pair, a 5-point 3-fold serial dilution standard curve was generated using cDNA template pools from six samples of the CT and WH treatment groups (n = 6, N = 12 total). No-template controls (NTC) were included to test for contamination and primer dimers. The 384-well format ViiA 7 Real-Time PCR system (Applied Biosystems/Life Technologies, Foster City, CA, USA) was used to perform amplification reactions in duplicate. qPCR amplifications were completed in 13 μL reactions with 1 × Power SYBR Green PCR Master Mix (Applied Biosystems/Life Technologies), 50 nM of each of the forward and reverse primers, and cDNA representing 10 ng of input total RNA as a starting point of the serial dilution. We applied the following real-time program: 1 cycle of 50°C for 2 min, 1 cycle of 95°C for 10 min, and 40 cycles of 95°C for 15 s and 60°C for 1 min, including fluorescence detection at the end of each 60°C step followed by a dissociation curve. Only primer pairs that had efficiencies between 84 and 108% (Additional file 11), and that generated an amplicon with a single melting curve without primer dimers were included for qPCR assays. All amplicons were examined through electrophoresis on 2% agarose gels along with a 1 Kb Plus DNA Ladder (TrackIt™, Invitrogen/Thermo Fisher, Carlsbad, CA, USA) to verify that the correct size fragment was amplified.

qPCR measurements (Fluidigm Biomark™)
The relative transcript expression values of 41 GOIs and the three previously evaluated normalizer genes were assessed for eight individual samples from each treatment group (n = 8, N = 24 in total) using the real-time qPCR Fluidigm Biomark™ HD system (Fluidigm, South San Francisco, CA, USA) based on 96.96 Dynamic Array™ IFC (GE-arrays) according to a protocol developed by Beemelmanns and Roth [144]. Firstly, a preamplification step was conducted for each sample by mixing 0.5 μL of a 500 nM STA primer pool (50 μM primer pair mix) with 2.5 μL of TaqMan-PreAmp Mastermix (Applied Biosystems, Waltham, MA, USA), 0.7 μL of nuclease-free water and 1.3 μL of cDNA (representing 200 ng of input total RNA), and using the following thermo-cycle protocol: 10 min at 95°C; 14 cycles of (15 s at 95°C, 4 min at 60°C). Then, the obtained PCR amplicons were diluted 1:10 with low EDTA-TE buffer. For the sample mix preparation, we combined 3.3 μL of preamplified cDNA with 3.5 μL of 2 × SsoFast EvaGreen Supermix with Low ROX (Bio-Rad Laboratories, Hercules, CA, USA) and 0.35 μL of 20 × DNA-Binding Dye Sample Loading Reagent (Fluidigm). The primer assay mix was prepared with 0.7 μL of 50 μM primer pair mix, 3.5 μL of Assay Loading Reagent (Fluidigm) and 3.15 μL of low EDTA-TE buffer. In a final step, 5 μL of each sample and assay mix was loaded on the 96.96 GEarray and fluorescence was measured using the Bio-mark™ HD system by running the GE Fast 96 × 96 PCR + Melt v2 thermal cycling protocol according to the manufacturer's instructions (Fluidigm). Transcript levels of the 44 GOIs were measured in two technical replicates, and we included two NTCs, two controls for genomic DNA contamination (no-reverse transcription), and two linker samples for inter-and intra-run calibrations.

qPCR data acquisition
For each technical replicate and sample, the mean threshold cycle (C T ), standard deviation (SD), and the coefficient of variation (CV) were calculated. As a quality control, C T values with a CV ratio greater than 4% [145] were removed from the dataset due to potential measurement errors. We performed geNorm analysis with the qBase+ software [143] based on the C T values of the three normalizer genes (eif3d, rpl32, pabpc1) [125,126,142] from all experimental samples. According to the geNorm analysis, the two normalizer genes rpl32 (geNorm M = 0.302) and eif3d (geNorm M = 0.313) were the most stable combination (geNorm V = 0.115). Based on the C T values obtained by the qPCR Fluidigm Biomark™ method, the two normalizer genes were stably expressed and showed mean C T differences below 0.438 ± 0.197 (mean ± SD) when comparing CT vs the WN or WH groups. Finally, the RQs of each gene were determined through normalization to the geometric mean (C T values) of rpl32 and eif3d expression, including the amplification efficiencies (Additional file 10), and setting the sample with the lowest normalized expression level as the calibrator sample (RQ value = 1.0) [143]. The corresponding FC values were calculated for each GOI by dividing its RQs by the mean of the control group.

Statistical analyses Principal component analysis (PCA)
All statistical tests were performed, and figures generated, in the R environment (v. 3.5.1) [146]. We applied multivariate statistical approaches to infer differences between the transcriptomic expression (log 2 ratios of microarray) of fish from the CT, WN and WH treatment groups based on: i) 17,072 detected probes; ii) 1111 DEPs shared between the WH and WN groups; iii) 789 DEPs specific for the WN group; and vi) 984 DEPs specific for the WH group. In addition, we assessed differential transcript expression patterns of qPCR genes (RQ values) grouped in broader categories of: i) 24 stress-related target genes (Themes # 1-4); and ii) 14 immune-related target genes (Theme # 5) ( Table 2). To explore differential gene expression patterns based on all of the abovementioned categories, we performed PCAs for graphical visualization using the dudi.pca function of the ade4 package in R [147]. For each PCA, the first two principal components (PC-1, PC-2) were plotted to obtain a projection of the whole dataset onto a small dimension and to account for the most relevant variance [148]. Then, the scores of PC-1 and PC-2 were extracted, and we fitted linear mixed-effect models for both PCs by applying the lme function implemented in the nmle package in R [149]. Statistical models were computed with the fixed interaction term 'treatment' and the random term 'tank' to account for between-tank variation (i.e., tank effect). Each model fit was graphically examined (histogram, qqplots) and residuals tested for normality (Shapiro-Wilk, p < 0.05). Finally, significant models were followed by Least-squares means post-hoc tests with Tukey's pvalue correction for multiple comparisons by applying the lsmeans function in R [150].

Gene-by-gene analysis
Differences in transcript expression between treatment groups were determined for each gene by fitting linear mixed-effect models and least-squares means post-hoc tests as explained above. Each model was graphically examined (hisImproved mitochondrial function in salmon togram, qqplots), model residuals were tested for normality (Shapiro-Wilk, p < 0.05), and RQ values were log 2 transformed. Seven outliers were removed (Bonferroni Outlier Test, p < 0.05) to fulfill assumptions.

Correlation analyses
First, component maps with the epPCA.inference.battery command of the package InPosition in R were computed [151] to: i) identify target genes of the highest importance; ii) relate the responses between genes; and to iii) assess the relationship between the expression of all 41 target genes and seven phenotypic traits (length, weight, CF, SGR, HSI, SSI and RVM). The incorporated battery of inference permutation tests calculated the component scores for each variable obtained by bootstrap ratios that are visualized in component maps [151]. To validate the gene expression results of the microarray and qPCR analyses, we correlated the mean log 2 FC values of the 41 target genes using Pearson correlations. Finally, to identify significant correlations between the expression (RQ) of the 41 genes and the seven phenotypic traits, we calculated a Pearson correlation matrix using the corrplot function in R [152].