The Transcriptomic Responses of Atlantic Salmon (Salmo salar) to High Temperature Stress Alone, and in Combination with Moderate Hypoxia

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

Download PDF

Research article

The Transcriptomic Responses of Atlantic Salmon (Salmo salar) to High Temperature Stress Alone, and in Combination with Moderate Hypoxia

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

This work is licensed under a CC BY 4.0 License

Journal Publication

published 12 Apr, 2021

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background: Increases in seawater 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 salmon’s hepatic transcriptome expression compared to control fish (CT: 12°C, normoxic) using 44K microarrays and qPCR.

Results: Overall, we identified 2,894 differentially expressed probes (DEPs, FDR < 5%), that included 1,111 shared DEPs, while 789 and 994 DEPs were specific to WN and WH fish, respectively. Pathway analysis suggested that the cellular mechanisms affected by the two experimental conditions were quite similar, with up-regulated genes functionally associated with 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 an impaired physiological and growth performance.

Conclusion: Increasing temperature to 20°C alone, and in combination with hypoxia, resulted in the up- and down-regulation of genes involved in similar important pathways in Atlantic salmon. However, the heat shock and immune responses of fish exposed to 20°C and hypoxia were more affected, and their transcriptional dysregulation was related to reduced performance. 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.

Epigenetics & Genomics

Climate Change

Increasing Temperature

Hypoxia

Transcriptomics

Biomarker Genes

Aquaculture

Temperature and oxygen are key environmental factors that influence the physiology, metabolism and survival of marine organisms, especially of ectothermic fish [1–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–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) in coastal regions [10, 11].

Thermal stress responses in fish have been widely investigated, and involve the expression of evolutionarily conserved stress-responsive 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–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–17, 20–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–4, 24]. Under hypoxic conditions, fish suppress energy-requiring processes like protein synthesis, aerobic metabolism, and mitochondrial energy production [25–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–27, 29–31]. To date, studies of fish physiology have reported large changes in specific aspects of the molecular responses to either hypoxia or increased temperature. Further, the fish in these experiments were generally exposed to an acute temperature increase and/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 ecosystems in temperate zones [33, 34]. Given the predicted increase in these two environmental stressors with global climate change [7, 9–11], it is of great importance that we assess the combined effects of these environmental challenges in experiments that simulate real-world 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–14 °C [37, 38]. In addition, the dissolved oxygen levels within the cages fluctuate substantially due to temperature, fish density, feeding and low water exchange [39–41], and frequent hypoxic events (~ 60–70% air saturation) occur in late summer [33, 34, 36]. These suboptimal conditions may negatively affect their 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]. Thus, it is critical that we: understand how processes at the molecular level are involved in regulating the physiological stress response and other important pathways; and identify putative biomarkers that can be used as tools to assess fish health and welfare, and can be incorporated into breeding programs with the goal of producing more robust fish.

In this study, we compared the hepatic transcriptional response of post-smolt salmon exposed to an incremental increase in temperature (12→20 °C at 1 °C week^− 1) alone or in combination with moderate hypoxia (~ 70% air saturation). The latter 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 that it has previously been shown to be an excellent target organ to characterize temperature and hypoxia stress response mechanisms in this fish species [14, 16]. An Agilent 44K salmonid oligonucleotide microarray platform [44] was employed in the current study 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 qPCR (Fluidigm Biomark™) on 41 target genes to: i) validate the microarray results; ii) examine the 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.

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. Further, all sections of this study adhere to the ARRIVE Guidelines for reporting animal research [45] and a completed ARRIVE guidelines checklist is included as Additional file 1.

Animal husbandry / initial holding

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, Newfoundland, Canada. Post-smolt Atlantic salmon (~ 1.5 years old) of Saint John River Strain origin (NB, Canada) obtained from Northern Harvest Sea Farms Ltd. were randomly distributed into six 2.2 m³ circular indoor fiberglass tanks receiving seawater (32 ppt salinity) at 15 L min^-1 (flow-through system). The fish were acclimated for four weeks under optimal conditions (~ 100–110% air saturation, 12 °C, 32 ppt salinity, 14 h light: 10 h dark photoperiod) and were 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). All fish in this experiment were implanted with Passive Integrated Transponder (PIT) tags (Loligo® Systems ISO 11784 certified PIT tags, Denmark) in their peritoneal cavity approximately two months before the experiment for identification purposes.

Experimental protocol

Three hundred and sixty Atlantic salmon with an initial average mass of 137.6 ± 1.3 g (mean ± S.E.) were randomly assigned to the following three treatment groups (two tank replicates, 60 fish tank^-1) as shown in Fig. 1: (1) Control (CT), constant temperature of 12 °C and ~ 100–110% air saturation for the duration of the experiment; (2) Warm & Normoxic (WN), incremental temperature increase (12→20 °C at 1 °C week^-1) at ~ 100–110% air saturation; and (3) Warm & Hypoxic (WH), decrease in oxygen content to ~ 70% air saturation (daily range ~ 65–75%) over one week, followed by two weeks of acclimation to this oxygen level, and then the same temperature regimen (12→20 °C at 1 °C week^-1) at ~ 70% air saturation. 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 daily monitored (YSI, ProODO, OH, USA) and ammonia and nitrite levels in the tanks were measured weekly (LaMotte test kit, Maryland, USA). During the experiment, salmon were carefully fed by hand to satiation twice daily (at 9:00 and 15:00) with the commercial salmon feed as indicated above. See Gamperl et al. [42] for more specific information on these experimental protocols.

For the current study we sampled eight fish 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). This minimum number of 24 samples was considered as sufficient to achieve statistical robustness and power (of 80%) to detect a significant effect (p < 0.05) with an estimated medium-large effect size (f² = 0.43), based on a priori multiple regression power calculation using the pwr.f2.test function of the pwr package in the R-software [46]. Further, this sample size was based on prior experiments with similar designs and treatments [42, 47].

For sampling, each fish was netted individually out of the tanks and euthanized in an aerated seawater bath (~ 10 L) containing a lethal dose (0.4 g L^-1) of the central nervous system depressant MS-222 (tricaine methanesulphonate; Syndel Laboratories, Nanaimo, BC, Canada) followed by a cerebral percussion (blow to the head between the eyes). We dissected 200–300 mg of liver tissue per fish, that was subsequently flash-frozen in liquid nitrogen and stored at -80 °C. Further, we measured the fish’s weight [g], fork length [cm], liver mass [g], spleen mass [g], and ventricle mass [g]. The growth and physiological indices of condition factor (CF), specific growth rate (SGR), hepato-somatic index (HSI), spleen-somatic index (SSI) and relative ventricular mass (RVM) were calculated, and these values (which are reported in Gamperl et al. [42]) were used in the current study for correlation analyses with gene expression.

RNA extraction, DNase treatment, and column purification

For RNA extraction, 100 mg of flash-frozen 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 system 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. [48] with minor changes. Crude RNA samples (100 µg, total amount) were mixed with an equal volume of Acid-Phenol:Chloroform (5:1 solution, pH 4.5, Ambion/Life Technologies, USA) and centrifugation was carried out at 17,000 g and 7 °C for 20 min. Then, 190 µL of the aqueous phase was precipitated with 0.1 volume of 3M 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, followed by centrifugation at 17,000 g and 7 °C for 30 min. The RNA pellets were then washed in 5 × volumes of 70% ethanol and centrifuged at 17,000 g and 7 °C for 20 min, air-dried 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, QIAGEN, Mississauga, ON, Canada) and 1 × of the manufacturer’s buffer. DNase-treated RNA samples were 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 for each of the treatment groups (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 of our group [49–52], and a sample-size and power calculation for microarrays using the sizepower function of the Bioconductor package in the R-software [53]. When considering a multiple-treatment design involving three treatment groups and the sample size of n = 6, we estimated an acceptable power of 88% to detect 2,000 differentially expressed genes.

An Agilent 44K Atlantic salmon oligonucleotide array platform (GEO accession: 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. [51] according to the MIAME guidelines [54]. 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. The quantity and quality of individual aRNA samples were evaluated using NanoDrop spectrophotometry and agarose gel electrophoresis, respectively. Thereafter, a common reference pool was generated by combining aRNA from all 18 samples (10 µg from each sample). Then, 20 µg of aRNA per individual sample or reference was precipitated and re-suspended in coupling buffer. The resulting aRNA was labeled with either Cy5 (individual samples) or with Cy3 (common reference) (GE HealthCare, Mississauga, ON, Canada) through a dye-coupling reaction, following the manufacturer’s instructions. The labeling efficiency and concentration of aRNA were determined using spectrophotometry (microarray function of NanoDrop), and the efficiencies ranged between 40–60 dye molecules per 1000 nt for all samples. Thereafter, 825 ng of labeled aRNA per sample was mixed with an equal amount of labeled common reference aRNA, and the resulting pools were fragmented and co-hybridized to 44K 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.

Microarray data acquisition

Each 4 × 44K Agilent microarray slide was scanned at 5 µm resolution with 90% of laser power using a ScanArray Gx Plus scanner and ScanExpress v4.0 software (Perkin Elmer, Waltham, MS, USA). 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₂-transformation and Loess normalization) were performed in R using mArray from the Bioconductor package [55] and with scripts adapted from Booman et al. [56] and Xue at al. [51]. 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 [57].

The whole microarray dataset consisted of 17,072 detected probes (normalized log₂ ratios, Cy5/C3 ratios) and is accessible on-line through the Gene Expression Omnibus (GEO) accession GSE146471 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146471).

Microarray data analysis

To identify significantly up- or down-regulated genes in response to the WN or WH treatments as compared to the control (CT) treatment, we applied permutational Significance Analysis of Microarrays (SAM) with a false discovery rate (FDR) of < 5% [58] using the implemented siggenes function of the Bioconductor package in R [59]. The SAM-identified differentially expressed probes (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 NCBI nr/nt database with E-value < 1 × 10^− 5 as filter criteria as described in Umasuthan et al. [49]. 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. [49].

A hierarchically clustered heat map was generated with normalized log₂ ratios for 2,894 SAM-identified DEPs (FDR < 5%) using the pheatmap R package [60] and applying Pearson correlation 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 differentially expressed genes (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) [61]. However, only non-redundant and annotated DEGs were recognized and included (WN group: 377 up-regulated DEGs and 798 down-regulated 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 two-sided hypergeometric test and the Benjamini-Hochberg p-value correction method [62]. The genes were mapped to Gene Ontology (GO) databases for biological, cellular, molecular and immune processes [63] using the Kyoto Encyclopedia of Genes and Genomes (KEGG) [64] and Reactome pathway databases [65] 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., down-regulated DEGs). The enriched leading 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 according 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 leading GO-terms of biological processes were summarised using the REVIGO cluster algorithm [66] and visualized with dot plots.

qPCR protocol

Gene selection and primer design

The selection of 41 microarray-identified genes of interest (GOIs) for validation purposes included significantly up- and down-regulated DEGs that were associated with important functional themes determined by previous GO/pathway term network analysis: Heat Shock Response (^#1), Cellular Stress (^#2), Oxidative Stress (^#3), Apoptosis (^#4), Immune Response (^#5), Transcription (^#7), and Cellular Metabolic Processes (^#10) (Table 1). We designed paralog-specific qPCR primers according to specifications detailed in Caballero-Solares et al. [50] 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 2. To identify any existing paralogs, and to verify the identity of each GOI, we performed BLASTn searches against the non-redundant nucleotide (nr/nt) and the expressed sequence tag (EST) databases of NCBI (year: 2018) (Additional file 3). For nine GOIs (camp-a, cat, cyp1a1, epx, gck, hif1α, hsp70, il8, and mmp9) the paralog-specific primer sequences were obtained from previous studies [67, 68] or the Genomic Applications Partnership Program (GAPP ^#6604) database, given that they were targeting the identical probe sequences (Additional files 2 and 3).

Table 1

Functional annotation and fold change (FC) transcript expression for 41 qPCR genes.
Gene Symbol¹	Gene Name²		Functional Theme³	Key Function (UniProt)⁴	FC WN⁵	FC WH⁶
Genes Related to Heat Shock Response
hsp70	Heat shock protein HSP 70		Heat Shock Response (^#1)	MF: Chaperone, Receptor; BP: Stress response	2.0 ± 0.37	2.0 ± 0.38
hsp90aa1	Heat shock protein HSP 90-alpha		Heat Shock Response (^#1)	MF: Chaperone; BP: Stress response	3.1 ± 0.31	3.7 ± 0.67
hsp90ab1	Heat shock protein HSP 90-beta		Heat Shock Response (^#1)	MF: Chaperone; BP: Stress response	1.3 ± 0.19	1.4 ± 0.22
hspd1	60 kDa heat shock protein		Heat Shock Response (^#1)	MF: Chaperone, Isomerase; BP: Stress response	0.3 ± 0.02	0.4 ± 0.04
serpinh1 (hsp47)	Serpin H1		Heat Shock Response (^#1)	MF: Chaperone; BP: Stress response, Collagen biosynthesis	5.5 ± 0.89	6.3 ± 0.78
Genes Related to Stress Response
calm (cam)	Calmodulin	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: Calcium ion binding; BP: Calcium-mediated signaling	0.4 ± 0.06	0.4 ± 0.04
cat	Catalase	Cellular Stress (^#2), Oxidative Stress (^#3) Apoptosis (^#4)		MF: Oxidoreductase, Peroxidase; BP: Hydrogen peroxide	1 ± 0.14	1 ± 0.11
cirbp	Cold-inducible RNA-binding protein	Heat Shock Response (^#1) Cellular Stress (^#2), Oxidative Stress (^#3)		MF: Activator, Repressor, RNA-binding; BP: Stress response	0.4 ± 0.03	0.3 ± 0.03
cldn3	Claudin 3	Cellular Stress (^#2),		MF: Transmembrane signaling; BP: Tight junction, Protein binding	0.8 ± 0.07	1 ± 0.1
cul3	Cullin 3	Cellular Stress (^#2), Apoptosis (^#4)		MF: Protein ubiquitination; BP: ER-Golgi transport, Transport	1.2 ± 0.08	1 ± 0.06
cyp1a1	Cytochrome P450 1A1	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: Oxidoreductase; BP: Oxidation-reduction process	0.2 ± 0.05	0.2 ± 0.04
egln2 (phd1)	Egl nine homolog 2	Cellular Stress (^#2), Oxidative Stress (^#3), Apoptosis (^#4)		MF: Oxidoreductase, Oxygen sensor activity; BP: Cell redox homeostasis, Response to hypoxia	0.4 ± 0.03	0.4 ± 0.05
gstt1	Glutathione S-transferase theta-1	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: Transferase; BP: Glutathione metabolic process	0.6 ± 0.09	0.5 ± 0.09
hcn1	Hyperpolarization-activated cyclic nucleotide-gated channel 1	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: Ion channel; BP: Ion transport, Transport	5.4 ± 1.27	5.7 ± 1.06
hif1α	Hypoxia-inducible factor 1 alpha	Cellular Stress (^#2), Oxidative Stress (^#3), Immune (^#5)		MF: DNA-binding; Oxygen sensor activity; BP: Transcription regulation, Response to hypoxia	0.60 ± 0.06	0.57 ± 0.06
igfbp2b1	Insulin-like growth factor-binding protein 2 - paralog b1	Oxidative Stress (^#3)		MF: Growth factor binding; BP: Growth regulation	0.6 ± 0.14	0.6 ± 0.15
jak2	Tyrosine-protein kinase JAK2-like	Cellular Stress (^#2), Oxidative Stress (^#3), Apoptosis (#4), Immune (^#5)		MF: Chromatin regulator, Tyrosine-protein kinase; BP: Apoptosis, Immunity	1.7 ± 0.22	1.5 ± 0.16
jund	Transcription factor Jun-D-like	Cellular Stress (^#2), Apoptosis (^#4), Immune (^#5)		MF: Activator, DNA-binding; BP: Transcription, Transcription regulation	3.2 ± 0.53	3 ± 0.35
ndufa1	NADH dehydrogenase 1 alpha subcomplex subunit 1	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: NADH dehydrogenase; BP: Electron transport, Respiratory chain	1 ± 0.05	1 ± 0.02
ndufa4	Cytochrome c oxidase subunit NDUFA4	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: NADH dehydrogenase; BP: Electron transport, Respiratory chain	2.2 ± 0.43	1.3 ± 0.45
prdx6	Peroxiredoxin 6	Cellular Stress (^#2), Oxidative Stress (^#3)		MF: Antioxidant, Peroxidase; BP: Cellular response to oxidative stress	0.4 ± 0.01	0.4 ± 0.02
rraga	Ras-related GTP-binding protein A	Cellular Stress (^#2), Oxidative Stress (^#3), Apoptosis (^#4)		MF: GTP-binding, Hydrolase; BP: Apoptosis	0.5 ± 0.04	0.5 ± 0.03
txn	Thioredoxin	Cellular Stress (^#2), Oxidative Stress (^#3) Apoptosis (^#4)		MF: Activator; BP: Oxidation-reduction process	1.9 ± 0.37	1.5 ± 0.26
ucp2	Mitochondrial uncoupling protein 2	Cellular Stress (^#2), Oxidative Stress (^#3), Apoptosis (^#4)		MF: Oxidative phosphorylation; BP: Transport	0.2 ± 0.04	0.2 ± 0.03
Genes Related to Immune Response
apod	Apolipoprotein D-like		Immune (^#5), Apoptosis (^#4), Oxidative Stress (^#3)	MF: Lipid-binding, Lipid transport; BP: Transport, Immunity	4.4 ± 2.73	7.7 ± 2.55
c1ql2	Complement C1q-like protein 2		Immune (^#5), Apoptosis (^#4)	MF: Protein binding; BP: Defence response to bacterium	6.4 ± 0.78	12.3 ± 4.47
c3	Complement C3-like		Immune (^#5)	MF: Signaling receptor binding; BP: Complement alternate pathway, Inflammatory response, Innate immunity	0.6 ± 0.11	0.6 ± 0.12
camp-a	Cathelicidin - paralog a		Immune (^#5)	MF: Antimicrobial peptide; BP: Bacterial defence	0.9 ± 0.2	3 ± 0.89
casp8	Caspase 8		Immune (^#5), Apoptosis (^#4)	MF: Hydrolase, Protease, Thiol-protease; BP: Apoptosis, Host-virus interaction	1.7 ± 0.12	1.6 ± 0.11
ctsh	Cathepsin H precursor		Immune (^#5), Apoptosis (^#4)	MF: Hydrolase, Protease, Thiol-protease; BP: Proteolysis	0.3 ± 0.03	0.3 ± 0.03
epx	Eosinophil peroxidase-like		Immune (^#5), Oxidative Stress (^#3)	MF: Oxidoreductase, Peroxidase; BP: Bacterial defence, Oxidative stress response	3.1 ± 0.85	2.3 ± 0.57
il8 (cxcl8)	Interleukin 8 (Chemokine CXC)		Immune (^#5)	MF: Cytokine; BP: Chemotaxis, Inflammatory response	1.4 ± 0.33	2.6 ± 0.73
irf2	Interferon regulatory factor 2		Immune (^#5)	MF: Activator, DNA-binding, Repressor; BP: Transcription regulation	0.8 ± 0.15	0.8 ± 0.17
mhcii (hla-dra)	MHC class ii antigen alpha chain		Immune (^#5)	MF: Peptide antigen binding; BP: Adaptive immunity, Immunity	1 ± 0.16	1.5 ± 0.32
mmp9	Matrix metalloproteinase 9		Immune (^#5), Apoptosis (^#4)	MF: Hydrolase, Metalloprotease, BP: Collagen degradation	2 ± 0.54	1.4 ± 0.25
nckap1l	Nck-associated protein 1-like		Immune (^#5), Apoptosis (^#4)	MF: Protein kinase activator; BP: Apoptosis	0.6 ± 0.05	0.6 ± 0.04
tapbp	Tapasin		Immune (^#5)	MF: Peptide antigen binding; BP: Immune response, Antigen presentation	0.8 ± 0.17	0.9 ± 0.21
tnfrsf6b	Tumor necrosis factor receptor superfamily member 6b		Immune (^#5), Apoptosis (^#4)	MF: Signaling Receptor; BP: Apoptosis	1.6 ± 0.45	2.8 ± 0.71
Genes Related to Cellular Metabolism
gck	Glucokinase		Cellular Metabolic Processes (^#10)	MF: Allosteric enzyme, Kinase, Transferase; BP: Glycolysis	0.5 ± 0.25	0.7 ± 0.35
pdk3	Pyruvate dehydrogenase [lipoamide] kinase isozyme 3		Cellular Metabolic Processes (^#10)	MF: Kinase, Transferase; BP: Carbohydrate metabolism, Glucose metabolism	1.8 ± 0.15	1.3 ± 0.11
Gene Related to Transcriptional Regulation
dnmt1	DNA (cytosine-5)- methyltransferase 1		Transcription (^#7)	MF: Chromatin regulator, DNA-binding, Methyltransferase; BP: DNA methylation	0.5 ± 0.02	0.5 ± 0.04
¹ Refers to the identity of the 44 K microarray-identified genes selected for qPCR validation. Gene abbreviations are
according to UniProt terminology. Further details about primer sequences and BLASTn hits for gene identification
purposes are given in Additional files 2 and 3.
² Categories of functional themes as identified in the GO/pathway term network analysis: 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) (see Figs. 3 and 4).
³ Key molecular function (MF) and biological process (BP) terms according to UniProt database for each gene.
⁴ Mean fold change ± S.E.M for the Warm & Normoxic (WN)group based on Fluidigm Biomark™ qPCR validation.
⁵ Mean fold change ± S.E.M for the Warm & Hypoxic (WH) group based on Fluidigm Biomark™ qPCR validation.

cDNA synthesis and primer efficiencies

First-strand cDNA templates for qPCR were synthesized in 20 µL reactions from 1 µg of DNaseI-treated, column-purified, total RNA with the QuantiTect^®Reverse - Transcription Kit (QIAGEN, Mississauga, ON, Canada) following the manufacturer’s protocol. The qPCR primer quality testing procedure was conducted according to previous protocols established by our group [50, 51, 56, 69]. 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 the WH treatment groups (n = 6, N = 12 total). No-template controls (NTC) were included to test for contamination and primer dimers. The ViiA 7 Real-Time PCR system with a 384-well format (Applied Biosystems/Life Technologies, Foster City, CA, USA) was utilized to perform amplification reactions in duplicate. Each of the PCR 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 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 attained standard curves with slopes of log₂ quantity vs threshold cycle (C_T) between − 3.8 and − 3.1, with efficiencies between 84–108% (Additional file 3), 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.

Selection of normalizer genes

Six normalizer genes [60S ribosomal protein L32 (rpl32, BT043656), eukaryotic translation initiation factor 3 subunit D (eif3d, GE777139), elongation factor 1-alpha 2 (ef1a12, BT058669), polyadenylate-binding protein 1 (pabpc1, EG908498), tubulin beta-2C chain (tubb2c, XM_014212152) and ATP binding cassette sub-family F member 2 (abcf2, BT071904)] from previous salmon transcriptome studies were selected as candidates [50, 51, 69]. The fluorescence threshold cycles (C_T) of eight samples per treatment group (n = 8, N = 24 in total) were measured for six normalizer genes following the method described in the previous paragraph with cDNA representing 5 ng of input total RNA. To identify the most stable normalizer, we conducted a geNorm analysis with the qBase + software [70] on the C_T values obtained from the 24 samples. Based on this analysis, the normalizer genes eif3d (geNorm M = 0.230), pabpc1 (geNorm M = 0.232) and rpl32 (geNorm M = 0.250) were identified as a stable combination (geNorm V = 0.108), and thus, included in the following qPCR measurements for all experimental samples.

qPCR measurements (Fluidigm Biomark^™)

The relative transcript expression values of 41 GOIs and the three previously evaluated reference genes were assessed for eight individual samples per three 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 [71]. Firstly, a pre-amplification step was conducted per 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 following thermo-cycle protocol: 10 min at 95 °C, 14 cycles: 15 s at 95 °C, 4 min at 60 °C. Then, the obtained PCR amplicons were 1:10 diluted with low EDTA-TE buffer. For the sample mix preparation, we combined 3.3 µL of pre-amplified 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 96.96 GE-array and fluorescence was measured using the Biomark^™ HD system by running the GE Fast 96 × 96 PCR + Melt v2 thermal cycling protocol according to the manufacturer’s instructions (Fluidigm). The transcript levels of 44 GOIs were measured in two technical replicates, while we included two NTCs, two controls for genomic DNA contamination (no-reverse transcription ‘no-RT’), and two linker samples for inter-run and between-run calibration. The 24 samples (Fish-ID ^#73 to ^#96) were assessed within a larger study (MICCSA) involving additional GOIs and RNA templates (N = 120 total), and the obtained threshold cycle (C_T) values for all samples of the MICCSA experiment are accessible on-line at the PANGAEA server (https://doi.org/10.1594/PANGAEA.913696).

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% [72] were removed from the dataset due to potential measurement errors. We performed geNorm analysis with the qBase + software [70] based on the C_T values of the three reference genes (eif3d, rpl32, pabpc1) [50, 51, 69] 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 considered to be the most stable combination (geNorm V = 0.115). Based on the C_T values obtained by the qPCR Fluidigm Biomark^™ method, the two endogenous control 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, respectively. Finally, the relative quantity (RQ) of each gene was determined through normalization to the geometric mean (C_T values) of the two endogenous reference genes (rpl32 and eif3d), including the amplification efficiencies (Additional file 2), and setting the sample with the lowest expression level as the calibrator sample (RQ value = 1.0) [70]. The corresponding fold change (FC) ratios were calculated using the raw RQs, and the mean of the control group was applied as a reference for each GOI.

Statistical analyses

Principal Component Analysis (PCA)

All statistical tests were performed, and figures generated, in the R environment (v. 3.5.1) [73]. We applied multivariate statistical approaches to infer differences between the transcriptomic expression (log₂ ratios of microarray) of fish from the CT, WN and WH treatment groups based on: i) 17,072 detected probes; ii) 1,111 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 1). To explore differential gene expression patterns based on all of the above-mentioned categories, we performed Principal Component Analysis (PCA) for graphical visualization using the dudi.pca function of the ade4 package in R [74]. 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 [75]. 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 [76]. 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 p-value correction for multiple comparisons by applying the lsmeans function in R [77].

Gene-by-gene analysis

Differences in gene expression between treatment groups were determined for each gene individually by fitting Linear mixed-effect models and Least-squares means post-hoc tests as explained above. Each model was graphically examined (histogram, qqplots), model residuals were tested for normality (Shapiro–Wilk, p < 0.05), and RQ values were log₂ 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 [78] to: i) identify target genes of the highest importance; ii) relate the responses between genes; and to iii) assess the relationship between the gene 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 [78]. To validate the gene expression results of the microarray and qPCR Fluidigm method, we correlated the mean log₂ fold change (FC) ratios of the 41 target genes using Pearson correlation. Finally, for the investigation of significant correlations between gene expression of 41 genes and seven phenotypic traits, we calculated a Pearson correlation matrix with the corrplot function in R [79].

Significance Analysis of Microarrays (SAM)

Based on the 17,071 detected microarray probes, the pairwise comparisons computed within SAM recognized 1,900 significantly differentially expressed probes (DEPs) for WN challenged fish and 2,105 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 4. When comparing the ‘WH vs CT’ and ‘WN vs CT’ DEP lists, we found that 1,111 DEPs (38%) were overlapping (i.e., were in common), 789 DEPs (28%) were WN-specific, and 994 DEPs (34%) were WH-specific (Fig. 2A). In summary, we identified a unique set of 2,894 DEPs when considering WN and/or WH conditions together, out of which 1,267 DEPs were up-regulated (WH = 732, WN = 135, shared for WH and WN = 400) and 1,627 DEPs were down-regulated (WH = 262, WN = 654, shared for WH and WN = 711) (Fig. 2A; Additional file 4). A hierarchically clustered heat map of the 2,894 DEPs displayed a robust cluster containing all CT samples that grouped separately from a 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 PCA based on 17,071 detected microarray features displayed a similar global transcriptional expression 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; CT vs WN p = 2e-4, CT vs WN p = 0.001, WN vs WN p = 0.489; 16.6% variation; Fig. 2C and Table 2). Similarly, we found a comparable global transcript expression response for the shared 1,111 DEPs with an equal cluster distribution for the WN and WH treatment groups as compared to the CT group (PC-1: p = 1e-04; CT vs WN p < 0.0001, CT vs WN p = 2e-04, WH vs WN p = 0.841; 61.5% variation; Fig. 2D and Table 2). However, the 789 WN-specific DEPs displayed a stronger differential expression as compared to the WN treatment, and this was indicated by a clear cluster separation between all three treatment groups (PC-1: p = 2e-4; CT vs WN p = 6e-4, CT vs WN p = 2e-4, WH vs WN p = 0.012; 46.6% variation; Fig. 2E and Table 2). On the contrary, the 994 WH-specific DEPs were more prominently featured in WH fish, as evidenced by an extreme shift away from the CT and WN groups (PC-1: p = 1e-4; CT vs WN p = 0.002, CT vs WN p < 0.0001, WH vs WN p = 0.007; 42.1% variation; Fig. 2F and Table 2).

Table 2

Temperature and hypoxia treatment effects on the extracted scores of the first two Principal Components.
Linear Mixed Effect Models Per Principal Components (PC-1 and PC-2)¹
Microarray²	PC	Var %	Factors	Num DF	Den DF	F-value	p-value	CT vs WN	CT vs WH		WN vs WH
All Probes (17,072)	PC-1	16.6	Intercept	1	12	20.6	7e-04
	PC-1	16.6	Treatment	2	3	384.6	2e-04	2e-04	0.001		0.489
	PC-2	10.4	Intercept	1	12	0.4	0.531
	PC-2	10.4	Treatment	2	3	1.0	0.477
DEPs WN&WH (1,111)	PC-1	61.5	Intercept	1	12	19.1	9e-04
	PC-1	61.5	Treatment	2	3	687.6	1e-04	< 0.0001	2 e-04		0.841
	PC-2	4.6	Intercept	1	12	0.3	0.623
	PC-2	4.6	Treatment	2	3	1.0	0.470
DEPs WN (789)	PC-1	46.6	Intercept	1	12	819.3	< 0.0001
	PC-1	46.6	Treatment	2	3	506.0	2e-04	6e-04	2e-04		0.012
	PC-2	6.3	Intercept	1	12	1.8	0.209
	PC-2	6.3	Treatment	2	3	0.9	0.506
DEPs WH (994)	PC-1	42.1	Intercept	1	12	252.3	< 0.0001
	PC-1	42.1	Treatment	2	3	875.0	1e-04	0.002	< 0.0001		0.007
	PC-2	7.7	Intercept	1	12	0.6	0.441
	PC-2	7.7	Treatment	2	3	0.0	0.982
qPCR³	PC	Var %	Factors	Num DF	Den DF	F-value	p-value	CT vs WN	CT vs WH	WN vs WH
All Target Genes (41)	PC-1	39.4	Intercept	1	18	147.5	< 0.0001
	PC-1	39.4	Treatment	2	3	318.4	3e-04	3e-04	6e-04		0.747
	PC-2	18.2	Intercept	1	18	46.9	< 0.0001
	PC-2	18.2	Treatment	2	3	1.2	0.416
Stress-Related Genes (24)	PC-1	45.9	Intercept	1	18	713.9	< 0.0001
	PC-1	45.9	Treatment	2	3	488.3	2e-04	1e-04	2e-04		0.537
	PC-2	19.6	Intercept	1	18	16.9	7e-04
	PC-2	19.6	Treatment	2	3	1.2	0.423
Immune-Related Genes (14)	PC-1	30.7	Intercept	1	18	1.5	0.237
	PC-1	30.7	Treatment	2	3	32.8	0.009	0.011	0.012		0.405
	PC-2	23.6	Intercept	1	18	0.0	0.870
	PC-2	23.6	Treatment	2	3	1.2	0.417
¹ Linear mixed effect models (lmes) were performed for PC-1 and PC-2 individually to assess the effect of the fixed
factor ‘treatment’ and including the random term ‘tank’. The variance explained by each PC is indicated in
percentage (%) values. For the pairwise comparisons between the Control (CT), Warm & Normoxic (WN), and
Warm & Hypoxic (WH) treatment groups, lsmeans post-hoc tests with Tukey's multiple comparisons adjustment of
p-values were applied. Significant p-values are marked in bold letters (p < 0.05).
² Statistical approach based on normalized log₂ ratios of 17,072 detected microarray probes, 1,111 WN&WH-specific
differentially expressed probes (DEPs), 789 WN-specific DEPs and 994 WH-specific DEPs.
³ Statistical approach based on the relative expression values (RQ) of 41 target genes, 24 stress-related genes and 14
immune-related genes measured with qPCR (Fluidigm Biomark™).

GO/pathway term network analysis

Salmon exposed to either the WN or WH treatments had similar patterns for enriched GO-terms and pathways (KEGG and/or Reactome) for up- and down-regulated differentially expressed genes (DEGs) (Figs. 3 and 4; Additional files 5–9). For both treatment groups, the up-regulated 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 5A, 5B, 6 and 7). On the contrary, the down-regulated DEGs for WN and WH fish were associated with a higher significance to enriched GO/pathway terms (p = 1e-5 to p = 1e-27) that were related to Proteolysis (^#8), Catabolic Processes (^#9), Cellular Metabolic Processes (^#10), and Oxidative Stress (^#3), and formed complex and interconnected functional clusters (Figs. 3B and 4B, Additional files 5C, 5D, 8 and 9). Below we list the shared and dissimilar enriched GO/term pathways for WN and WH fish visualized in 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 6–9.

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 larger cluster (~ 8 terms) with ‘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) as the most significantly enriched leading terms (Figs. 3A; Additional file 6). On the contrary, WH fish had two larger and highly interconnected clusters (~ 30 terms) associated with the chaperone-mediated heat shock response and including enriched leading terms: ‘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-HSA:381119) (Fig. 4A; Additional file 7). Also, similar GO-terms were enriched in WN and WH fish related to ‘chaperone-mediated autophagy processes’ (WN: GO:0061684; WH: GO:0061741), ‘chaperone and protein folding responses’ (WN: GO:0061077; WH: GO:0030968), and ‘regulation of tau-protein kinase activity’ (WN + WH: GO:1902947, GO:1902949) (Additional file 5A-B).

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) represented by the following enriched leading 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) (Fig. 3A; Additional file 6). Whereas, for fish subjected to WH conditions, we identified the following unique oxidative stress-related leading terms: ‘antioxidant activity’ (GO:0016209), ‘positive regulation of release of cytochrome c from mitochondria’ (GO:0090200), ‘positive regulation of catecholamine metabolic process’ (GO:0045915), ‘positive regulation of transcription in response to endoplasmic reticulum stress’ (GO:199044), and ‘thioredoxin oxidation’ (R-HSA:73646) (Fig. 4A; Additional file 7).

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’ (GO:1990442), ‘regulation of nitrosative stress-induced intrinsic apoptotic signaling’ (GO:1905258), ‘modulation of programmed cell death / apoptotic process in other organism’ (GO:0044531, GO:0044532), and ‘regulation of hydrolase activity’ (GO:0051336) (Fig. 3A; Additional file 6). On the contrary, WH fish had five smaller clusters containing ~ 10 enriched GO/pathway terms associated with apoptotic processes like ‘regulation of intrinsic apoptotic signaling pathway’ (GO:2001242), ‘positive regulation of endothelial cell apoptotic process’ (GO:2000353), and ‘transcription factor AP-1 complex’ (GO:0035976) (Fig. 4A; Additional file 7). Also, similar GO-terms related to ‘cysteine-type endopeptidase activity involved in apoptotic signaling pathway’ (GO:2001267, GO:0097199) and ‘regulation of nitrosative stress-induced intrinsic apoptotic signaling pathway’ (GO:1905258, GO:1990442, GO:1905259) were enriched in WN and WH fish (Figs. 3A and 4A; Additional files 5A-B, 6 and 7).

Immune Response (^#5)

Pronounced differences between WH and WN groups were found for pathways related to immune response (^#5) (Figs. 3A and 4A; Additional files 6 and 7). Up-regulated DEGs in the WN group were related to enriched immune-related pathways (five groups, ~ 13 terms) such as ‘IL-17 signaling pathway’ (KEGG:04657), ‘Interleukin-4 and Interleukin-13 signaling’ (R-HSA:6785807), ‘monocyte chemotactic protein-1 production’ (GO:0071605), ‘activation of matrix metalloproteinases’ (R-HSA:1592389), and ‘regulation of blood coagulation’ (GO:0030193) (Fig. 3A; Additional file 6). On the contrary, the up-regulated DEGs of the WH group formed more complex immune-related functional clusters (12 groups, ~ 43 terms) that were associated with ‘neutrophil/granulocyte migration and granulation’ (GO:0097530, GO:0071621, GO:1990266, GO:0030593), ‘phagosome’ (KEGG:04145), ‘JAK2 growth hormone receptor signaling’ (R-HSA:6784189) and various interleukin signaling pathways such as ‘IL-17 signaling pathway’ (KEGG:04657) and ‘IL-4 and IL-13 signaling’ (R-HSA:6785807) (Fig. 4A; Additional files 5B and 7). GO/pathway terms connected to anti-viral responses were also activated such as ‘interferon signaling’ (R-HSA:913531), ‘herpes simplex virus 1 infection’ (KEGG:05168) and ‘MHC class I protein complex assembly’ (GO:0002397). Further, one cluster connected to extracellular matrix processing such as ‘collagen type III degradation by MMP1,8,9,13’ (R-HSA:1474213) and also ‘blood coagulation extrinsic pathways’ (GO:0007598) was enriched (Fig. 4A; Additional files 5B and 7).

Protein Processing & Localization (^#6) and Transcription (^#7)

Both experimental groups had complex functional network clusters linked to Protein Processing & Localization (^#6) and Transcription (^#7) (Figs. 3A and 4A) including ‘translation initiation’ (GO:0006413, R-HSA:72613), ‘protein localization to endoplasmic reticulum’ (WN: GO:0070972, WH: GO:0072599), ‘mRNA metabolic process’ (WN: GO:0016071, WH: GO:0000184) and ‘spliceosome’ or ‘RNA metabolic processes’ (WN: GO:0000184, WH: GO:0043484, KEGG:03040) (Figs. 3A and 4A; Additional files 5A-B, 6 and 7).

GO/pathway term networks associated with down-regulated DEGs

Proteolysis (^#8), Catabolic Processes (^#9), and Cellular Metabolic Processes (^#10)

The most highly enriched terms were related to ‘small molecule metabolic process’ (GO:0044281/ GO:0044283; WN: p = 1e-27; WH: p = 3e-10), ‘carboxylic acid metabolic process’ (GO:0019752; WN p = 1e-21, WH p = 1e-15), ‘organic acid metabolic process’ (GO:0006082; WN p = 2e-20, WH p = 3e-15), and ‘proteasome’ (KEGG:03050; WN p = 9e-20, WH p = 2e-13) (Figs. 3B and 4B; Additional files 8 and 9). For both treatment groups, we found large clusters with functional terms linked to Proteolysis (^#8) and Catabolic Processes (^#9) such as ‘proteasome’ (KEGG:03050, GO:0000502) and ‘peptidase complex’ (GO:1905368) (Figs. 3B and 4B; Additional files 8 and 9). The majority of enriched metabolic pathways can broadly be categorized as Cellular Metabolic Processes (^#10) [e.g., ‘cellular amino acid metabolic process’ (GO:0006520), ‘carboxylic acid metabolic process’ (GO:0019752), ‘lipid metabolic process’ (GO:0006629), ‘fatty acid metabolic process’ (GO:0006631), ‘nucleobase-containing small molecule metabolic process’ (GO:0055086), ‘glucose 6-phosphate metabolic process’ (GO:0051156), and ‘citrate cycle’ (TCA cycle) (KEGG:00020, GO:0006101) (Figs. 3B and 4B; Additional files 5C-D, 8 and 9). Interestingly, the down-regulated DEGs of WH fish formed four interconnected clusters associated with Oxidative Stress (^#3) like ‘oxidoreductase activity’ (GO:0016491), ‘oxidoreductase activity, acting on CH-OH group of donors’ (GO:0016614) and ‘mitochondrion’ (GO:0005739) (Fig. 3B; Additional file 9), while WN fish also had two smaller clusters connected to ‘oxidation-reduction processes’ (GO:0055114, GO:0016614) and ‘mitochondrion’ (GO:0005739) (Fig. 4B; Additional file 8).

qPCR validation of microarray approach

We chose a subset of 41 treatment-responsive probes to confirm the microarray results using a gene-targeted qPCR approach. A significant positive correlation (R = 0.87; p < 1.7e-13) between the log₂ FC means of the 41 genes measured via the 44K microarray and qPCR methods indicates that the microarray results were validated with great confidence (Additional file 10). 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 11). Minor differences were found for the genes camp-a, cat and cldn3, whilst igfbp2b1, irf2, nckap1l and tapbp had a difference in the direction of expression change (Additional file 11). Out of the 41 GOIs, we found 23 with a significant ‘treatment’ effect in the statistical model (Table 3) and a selection of GOIs are illustrated in Fig. 5. 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 cross-hybridization on the 44K array.

Table 3

Temperature and hypoxia treatment effects on the transcript expression of 41 target genes (qPCR validation).
Linear Mixed Effect Models Per Target Gene¹
	Intercept		Treatment		Lsmeans Post-Hoc Test
Gene	F-value	p-value	F-value	p-value	Lsmeans Post-Hoc Test
apod	165.5	< 0.0001	15.46	0.026	CT-WH
c1ql2	1954.6	< 0.0001	48.87	0.005	CT-WN, CT-WH
c3	132.1	< 0.0001	2.53	0.227
calm	1875.4	< 0.0001	84.92	0.002	CT-WN, CT-WH
camp-a	429.4	< 0.0001	6.02	0.089
casp8	249.8	< 0.0001	16.97	0.023	CT-WN, CT-WH
cat	19704.0	< 0.0001	4.54	0.124
cirbp	356.3	< 0.0001	84.34	0.002	CT-WN, CT-WH
cldn3	739.8	< 0.0001	13.75	0.031	WN-WH
ctsh	509.5	< 0.0001	76.82	0.003	CT-WN, CT-WH
cul3	1087.0	< 0.0001	2.28	0.250
cyp1a1	362.9	< 0.0001	30.48	0.010	CT-WN, CT-WH
dnmt1	1240.3	< 0.0001	40.65	0.007	CT-WN, CT-WH
egln2	5750.5	< 0.0001	76.76	0.003	CT-WN, CT-WH
epx	245.1	< 0.0001	19.68	0.019	CT-WN, (CT-WH)
gck	189.4	< 0.0001	2.47	0.232
gstt1	802.0	< 0.0001	18.43	0.021	(CT-WN), CT-WH
hcn1	225.1	< 0.0001	41.89	0.006	CT-WN, CT-WH
hif1α	394.9	< 0.0001	10.88	0.042	(CT-WN), CT-WH
hsp70	61.6	< 0.0001	3.18	0.181
hsp90aa1	291.6	< 0.0001	45.56	0.006	CT-WN, CT-WH
hsp90ab1	255.8	< 0.0001	2.05	0.275
hspd1	1795.4	< 0.0001	167.12	0.001	CT-WN, CT-WH
igfbp2b1	81.7	< 0.0001	1.68	0.323
il8	124.1	< 0.0001	3.72	0.154
irf2	403.4	< 0.0001	0.89	0.499
jak2	318.8	< 0.0001	8.01	0.053
jund	166.4	< 0.0001	16.66	0.024	CT-WN, CT-WH
mhcii	93.8	< 0.0001	1.69	0.322
mmp9	76.1	< 0.0001	9.23	0.052
nckap1l	126.7	< 0.0001	17.66	0.022	CT-WN, CT-WH
ndufa1	4303.5	< 0.0001	0.58	0.611
ndufa4	93.1	< 0.0001	8.26	0.060
pdk3	1163.6	< 0.0001	16.22	0.025	CT-WN, (CT-WH)
prdx6	1814.2	< 0.0001	158.74	0.001	CT-WN, CT-WH
rraga	1624.3	< 0.0001	46.90	0.006	CT-WN, CT-WH
serpinh1	54283.9	< 0.0001	31.67	0.010	CT-WN, CT-WH
tapbp	69.2	< 0.0001	0.14	0.878
tnfrsf6b	69.6	< 0.0001	4.61	0.122
txn	82.9	< 0.0001	3.21	0.180
ucp2	409.9	< 0.0001	147.53	0.001	CT-WN, CT-WH
	NumDF = 1, DenDF = 18		NumDF = 2, DenDF = 3
¹ Linear mixed effect models (lmes) were performed for each gene individually to assess the effect of the
fixed factor ‘treatment’ and including the random term ‘tank’. To investigate pairwise comparisons between the
Control (CT), Warm & Normoxic (WN) and Warm & Hypoxic (WH) treatment groups the significant models were
followed by lsmeans post-hoc test with Tukey's multiple comparisons adjustment of p-values. Significant values are
marked in bold letters (p < 0.05) and trends are indicated in italics and/or in brackets (0.05 < p < 0.1).
Additional files
Additional file 1: Checklist. Completed “ARRIVE Guidelines Checklist” for reporting animal data in this study. (PDF 451.9 KB)
Additional file 2: Primers used for qPCR validation analyses. Listed are gene symbol, gene name, the sequence of forward and reverse primers, the amplicon size in base pairs, GenBank accession number used for primer design, and the determined primer efficiencies. See Additional file 3 for further details about BLASTn hits and primer properties. (DOCX 34.0 KB)
Additional file 3: Details about target genes and primers used for qPCR validation analyses. Listed are specific details about the qPCR genes and primers, BLASTn results of amplicon sequences, and primer properties for 41 target genes and two normalizer genes. (XLSX 43.0 KB)
Additional file 4: Results of Significance Analysis of Microarrays (SAM, FDR < 5%). Listed are the fold change (FC) values for significantly differentially expressed probes (DEPs) that were identified in the liver of Atlantic salmon exposed to (A) high temperature and normoxia (WN: 20 °C, 100% air sat.) or (B) high temperature and hypoxia (WH: 20 °C, ~ 70% air saturation) as compared to control conditions (CT: 12 °C, 100% air sat.) (n = 6, N = 18 total). (XLSX 15 MB)
Additional file 5. Enrichment GO-term dot plot for up- and down-regulated differentially expressed genes (DEGs) in Atlantic salmon that were subjected to Warm & Normoxic (WN: 20 °C, 100% air sat.) or Warm & Hypoxic (WH: 20 °C, ~ 70% air sat.) conditions. The dot plots represent non-redundant significantly enriched gene ontologies (GO) terms of biological processes after application of REVIGO's redundancy elimination algorithm for (A) up-regulated DEGs of the Warm & Normoxic group; (B) up-regulated DEGs of the Warm & Hypoxic group; (C) down-regulated DEGs of the Warm & Normoxic group; and (D) down-regulated DEGs of the Warm & Hypoxic group. The colour scheme corresponds to the log₁₀ adjusted p-values (Benjamini and Hochberg method), and the diameter of the dots represents the number of DEGs that were identified to be significantly associated with this specific term. (PDF 271.8 KB)
Additional file 6: Results of GO/pathway term network analysis (ClueGO) associated with up-regulated differentially expressed genes (DEGs) of the Warm & Normoxic (WN) treatment group. Listed are the significantly enriched gene ontology (GO) terms, and KEGG or Reactome pathways (hypergeometric test p < 0.05 with Benjamini-Hochberg correction). GO/pathway term annotations were obtained using the GO database for biological process, cellular component, molecular function, and immune processes, KEGG and Reactome pathway databases. The enriched GO/pathway terms are ordered according to 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). The enriched GO/pathway terms are illustrated in corresponding network Fig. 3A. (XLSX 331.0 KB)
Additional file 7: Results of GO/pathway term network analysis (ClueGO) associated with up-regulated differentially expressed genes (DEGs) of the Warm & Hypoxic (WH) treatment group. Listed are the significantly enriched gene ontology (GO) terms, and KEGG or Reactome pathways (hypergeometric test p < 0.05 with Benjamini-Hochberg correction). GO/pathway term anotations were obtained using the GO database for biological process, cellular component, molecular function, and immune processes, KEGG and Reactome pathway databases. The enriched GO/pathway terms are ordered according to 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). The enriched GO/pathway terms are illustrated in corresponding network Fig. 4A. (XLSX 93.9 KB)
Additional file 8: Results of GO/pathway term network analysis (ClueGO) associated with down-regulated differentially expressed genes (DEGs) of the Warm & Normoxic (WN) treatment group. Listed are the significantly enriched gene ontology (GO) terms, and KEGG or Reactome pathways (hypergeometric test p < 0.05 with Benjamini-Hochberg correction). GO/pathway term annotations were obtained using the GO database for biological process, cellular component, molecular function, and immune processes, KEGG and Reactome pathway databases. The enriched GO/pathway terms are ordered according to 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). The enriched GO/pathway terms are illustrated in corresponding network Fig. 3B. (XLSX 91.8 KB)
Additional file 9: Results of GO/pathway term network analysis (ClueGO) associated with down-regulated differentially expressed genes (DEGs) of the Warm & Hypoxic (WH) treatment group. Listed are the significantly enriched gene ontology (GO) terms, and KEGG or Reactome pathways (hypergeometric test p < 0.05 with Benjamini-Hochberg correction). GO/pathway term annotations were obtained using the GO database for biological process, cellular component, molecular function, and immune processes, KEGG and Reactome pathway databases. The enriched GO/pathway terms are ordered according to 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). The enriched GO/pathway terms are illustrated in corresponding network Fig. 4B. (XLSX 89.5 KB)
Additional file 10: Pearson correlation between mean log₂ fold change (FC) ratios of microarray probes and log₂ fold changes (FC) of 41 target genes measured with qPCR (Fluidigm Biomark™). The relationship between data obtained for the 41 target genes by both methods was estimated using Pearson’s product-moment correlation. The correlation coefficient R and the significance of correlation (p < 0.0001***) are presented in the figure. (PDF 31.7 KB)
Additional file 11: Comparison of fold change (FC) values for 41 target genes obtained by 44K microarray (Agilent) and qPCR (Fluidigm Biomark™) approach. (DOCX 28.4 KB)
Additional file 12: Bootstrap ratios from the first two components (CP-1 and CP-2) of the component maps based on 41 target genes and seven phenotypic traits for the Warm & Normoxic (WN) and Warm & Hypoxic (WH) treatment groups. Bootstrap ratios were performed to identify the variables that contributed significantly to the variance of a given component. Bold values indicate bootstrap ratios whose magnitude exceed +/-2 and are considered as significant. (DOCX 22.0 KB)

Gene-by-gene analysis

Genes related to Heat Shock Response (^#1)

The heat shock response-related genes serpinh1 (p = 0.010; Fig. 5A) and hsp90aa1 (p = 0.006; Fig. 5B, Table 3) were significantly higher expressed in fish challenged with the WN and WH conditions in comparison to CT conditions. 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 up-regulated in WN and WH fish as compared to CT fish (p = 0.006; Fig. 5D, Table 3). On the contrary, 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 (p = 0.010; Fig. 5G), egln2 (p = 0.003; Fig. 5H), prdx6 (p = 0.001; Fig. 5I), rraga (p = 0.006; Fig. 5J) and ucp2 (p = 0.001; Fig. 5K) (Table 3). Two stress- and hypoxia-related genes, gstt1 (p = 0.021; Fig. 5L) and hif1α (p = 0.042; Fig. 5M), were only significantly down-regulated in the WH group as compared to the CT group (Table 3). The gene dnmt1 related to 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 connected to metabolic processes was significantly up-regulated 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)

Three genes related to apoptotic processes, casp8 (p = 0.023; Fig. 5P), jund (p = 0.024; Fig. 5Q) and jak2 (p = 0.053; Fig. 5R) were more highly expressed in fish challenged with the WN and WH conditions in comparison to the CT group (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)

We found higher expression levels in the treatment groups, as compared to the CT group, for three immune-related genes: apod (p = 0.026; Fig. 5T), c1ql2 (p = 0.005; Fig. 5U) and epx (p = 0.019; Fig. 5V). However, they had different expression profiles (Table 3). The expression of apod was only significantly higher in the WH group, while the expression of epx was only higher in the WN group as compared to the CT group (Table 3). In contrast, the gene c1ql2 was significantly up-regulated in both experimental groups as compared to the CT group (Table 3). The expression of nckap1l was significantly down-regulated in the WN and WH group as compared to the CT group (p = 0.022; Fig. 5W), this result was contrary to the up-regulated expression of the microarray probe (Additional file 11).

PCA analysis on 24 stress-related (^#1–4) and 14 immune-related (^#5) genes

The PCA analysis based on 24 stress-related genes revealed similar transcript expression profiles for salmon exposed to WN and WH conditions in comparison to the CT condition (PC-1: p = 2e-04; CT vs WN p = 1e-04, CT vs WN p = 2e-04, WH vs WN p = 0.537; 45.9% variance; Fig. 6A and Table 2). A similar pattern was found for the global dysregulation of 14 immune-related genes. However, the WH group had a wider dispersion of variance as compared to both the CT and WN groups (PC-1: p = 0.009; CT vs WN p = 0.011, CT vs WN p = 0.012, WH vs WN p = 0.405; 30.7% variance; Fig. 6B and Table 2).

Correlation between gene expression and phenotypic traits

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, 37.4%) (Fig. 6C). The following up-regulated genes correlated significantly with CP-1: c1ql2 (CP-1= -6.3), hsp90aa1 (CP-1= -5.3), casp8 (CP-1= -5.2), serpinh1 (CP-1= -5.0), pdk3 (CP-1= -4.3), jund (CP-1= -3.6), hcn1 (CP-1= -3.6), jak2 (CP-1= -3.1), hsp70 (CP-1= -2.6), ndufa4 (CP-1= -2.3) and epx (CP-1= -2.3) (Fig. 6C; Additional file 12). On the contrary, the down-regulated genes that showed a strong association with PC-1 were: prdx6 (CP-1 = 11.6), egln2 (CP-1 = 8.9), cirbp (CP-1 = 8.5), hspd1 (CP-1 = 8.5), ctsh (CP-1 = 8.3), dnmt1 (CP-1 = 7.8), rraga (CP-1 = 7.2), calm (CP-1 = 6.9), ucp2 (CP-1 = 5.7), cyp1a1 (CP-1 = 5.1), hif1α (CP-1 = 4.3) and gstt1 (CP-1 = 4.3) (Fig. 6C, Additional file 12). Since the main treatment effect was explained by PC-1 (Table 2), all the genes with significant associations to CP-1 are likely to be key drivers of the phenotypic changes in WN fish (reported in Gamperl et al. [42]). Out of the seven phenotypic characteristics for the WN group, only SGR correlated significantly with CP-2 (Fig. 6C; Additional file 12), and six genes showed positive associations (R between 0.5 and 0.65) with SGR (irf2, tapbp, hsp70, hsp90aa1, igfbp2b1, and ndufa1; Fig. 6E).

Correlation analysis for WH fish

In the component map based on the response variables obtained for WH fish (Fig. 6D), 19 up-regulated and 22 down-regulated genes were segregated on opposing directions along CP-1 (CP-1, 38.5%). For the WH group we found similar up-regulated genes as compared to WN group that correlated with the main explanatory CP-1: serpinh1 (CP-1=-6.8), casp8 (CP-1=-4.6), jund (CP-1= -4.5), hcn1 (CP-1= -4.0), hsp90aa1 (CP-1= -3.2), hsp70 (CP-1= -3.0), jak2 (CP-1= -2.9), c1ql2 (CP-1= -2.7), epx (CP-1= -2.3) and pdk3 (CP-1= -2.1) (Fig. 6D; Additional file 12). Interestingly, only in WH fish four immune response-related genes were significantly correlated with CP-1: apod (CP-1= -2.5), tnfrsf6b (CP-1= -2.3), il8 (CP-1= -2.1), and camp-a (CP-1= -2.0) (Fig. 6D; Additional file 12). We also identified the same set of down-regulated genes for WH fish, but with a high contribution: prdx6 (CP-1 = 10.3), rraga (CP-1 = 8.6), ctsh (CP-1 = 8.0), cirbp (CP-1 = 7.9), calm (CP-1 = 7.4), hspd1 (CP-1 = 7.0), egln2 (CP-1 = 7.0), ucp2 (CP-1 = 6.1), dnmt1 (CP-1 = 5.9), cyp1a1 (CP-1 = 5.2), gstt1 (CP-1 = 4.6) and hif1α (CP-1 = 4.1) (Fig. 6D; Additional file 12). Numerous associations were also evident between gene expression and growth indices for WH fish. For example weight (CP-1 = 3.1), CF (CP-1 = 3.1), SGR (CP-1 = 2.7), and length (CP-1 = 2.0) (Fig. 6D; Additional file 12) correlated positively with the expression of down-regulated genes. When using Pearson correlation coefficients, 12 down-regulated genes (ctsh, nckap1, cirbp1, hspd1, calm, egln2, hif1α, ucp2, gstt1, prdx6, rraga and dnmt1) showed strong positive correlations (R between 0.51 to 0.81) with reduced growth parameters (i.e., weight, length, CF or SGR) as previously published in Gamperl et al. [42] (Fig. 6F). On the other hand, seven up-regulated genes (apod, casp8, camp-a, il8, hsp90aa1, serpinh1 and jund) correlated negatively (R between − 0.52 to -0.73) with these phenotypic characteristics (Fig. 6F). Finally, the immune-related genes apod (R = 0.59) and camp-a (R = 0.57) correlated positively with SSI (Fig. 6F).

Seasonal fluctuations in water temperature and oxygen content are unavoidable changes 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 the occurrence of suboptimal temperatures (both high and low) and low water oxygen levels (hypoxia) may become more frequent and severe in coastal areas [7–11], and ultimately more challenging for aquaculture species [80]. 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.) that realistically simulated summer conditions in salmon aquaculture sea-cages [33, 34, 41]. Then, functional genomic approaches were employed (Agilent 44K salmonid oligonucleotide microarray and qPCR Fluidigm Biomark^™) 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 similar global gene expression responses that were highly distinctive from the control (12 °C and normoxic) fish. Overall, we identified a set of 2,894 DEPs, out of which 1,111 DEPs (38%) were shared between the WH and WN treatment groups. The 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 large 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 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], and this reinforces the biological relevance of these genes and pathways and 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 cells and consequently tissues from structural damage [81–83]. 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 the loss of functionality [82, 83]. In addition, HSPs have a wider role in relation to the function of the immune system (i.e., antigen presentation) [84, 85], apoptosis [86] and protection from oxidative stress [87]. HSPs are essential regulators of the cellular stress response in aquatic ectotherms [88], and their expression is well characterized in teleost species as they function as protective proteins against acute and chronic thermal stress [13–17, 21–23] and during hypoxia [26, 29, 30]. After the incremental temperature challenge to 20 °C (alone or combined 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, Serpin H1) was one of the most up-regulated target genes (qPCR validated) for both warmed groups in comparison to the control group, and this is in agreement with previous studies on over ten fish species [89]. Serpin H1 binds very specifically to collagens and procollagens to facilitate their assembly and stabilization and plays an important role in collagen biosynthesis [90]. Moreover, Serpin H1 is involved in the breakdown of reactive oxygen species (ROS) produced during oxygen stress as recently shown in rainbow trout (Oncorhynchus mykiss) [91]. 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 generated ROS to maintain cellular homeostasis during this thermal challenge. Likewise, the expression of the gene hsp90aa1 (alias hsp90-alpha, Heat Shock Protein 90 alpha) that codes for the inducible form of HSP90 and its constitutive counterpart hsp90ab1 (alias hsp90-beta, Heat Shock Protein 90 beta) were both up-regulated 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, 92]. 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 [88, 93, 94]. In our study, hsp90aa1 was interconnected with most of the HSP-related GO/pathway terms like ‘HSF1 activation’ and ‘dissociation of cytosolic HSF1:HSP90 complex’, 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, the expression of this transcript was positively correlated with hsp90aa and serpinh1 expression levels, and was associated with enriched GO/pathway terms such as ‘HSP70:HSP40:nascent protein’. 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 [88, 89]. Thus, in the current study, inducible forms of HSP70 may have been involved in molecular chaperone processes within the liver cells of stress exposed fish to assist with the folding of nascent polypeptide chains and the repair and degradation of altered or denatured proteins [95]. 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 point to the presence of ER-stress and UPR. The ER is considered to be a factory for secretory proteins with quality-control systems ensuring the correct folding of proteins and their vesicular cellular transport [96]. The ER-stress response initiates the UPR to prevent the assimilation of unfolded proteins and restore ER function [96].

In summary, we found a highly active cellular HSP response and chaperone-mediated ER-stress response 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 a larger and more interconnected cluster associated with the HSP response (Theme ^#1) and different unique pathways (e.g., ‘UPR’). These latter results indicate that there was a greater induction of these essential cell maintenance processes in the liver of WH fish during this climate change scenario.

Apoptosis and programmed cell death

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 [97]. High temperature and hypoxia can induce oxidative stress that mediates apoptotic processes in fish tissues [20, 98–101]. 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 with 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 pathways (e.g., unfolded protein binding) and immune pathways (e.g., MHC protein complex binding) that were linked with apoptotic processes. This suggests that enhanced HSP transcript expression (i.e., hspaa1 and hspd8) might have been associated with cell resistance to apoptotic cell death [86].

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 (Caspase 8) was the most frequently occurring transcript in enriched apoptosis-related GO/pathways 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 [102]. On the contrary, the significantly up-regulated transcript jund (Transcription factor JunD) in WN and WH fish indicates that processes were activated to protect cells from senescence or apoptosis by acting as a modulator of the p53 signaling cascade [103].

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 other pathways likely prevented extensive programmed cell death which could have resulted in hepatic necrosis or impacted liver function [1, 97, 100, 101]. 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. However, at present, we do not have a reasonable explanation for this finding.

Immune response

Several studies have shown that exposure to elevated temperatures [15, 17–19, 22] or hypoxia [6, 52] affects immune-related gene expression in fish. However, these studies have 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 systems (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 the 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 signaling’), and antiviral responses (i.e., ‘herpes simplex virus 1 infection’). In addition, the 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 dysregulation of several genes only in WH fish such as apod, camp-a, c1ql2, il8, mhcii, and tnfrsf6b, and some of these correlated strongly with health parameters. For instance, apod (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 showed a positive correlation with spleen-somatic index, highlighting its potential relevance for fish health. The extracellular glycoprotein APOD has multiple functions that involve the immune response, chemoreception, proteolysis and lipid oxidation [104], and potentially plays a fundamental role in cellular-stress and immune responses. The higher expression of c1ql2 (C1q-like Protein 2) transcripts (WN: 6.4-fold; WH: 12.3-fold, qPCR validated) suggests that the classical complement system pathway was activated. This is an essential innate defence mechanism in teleost fish that detects and destroys invading pathogens by bacteriolysis [105]. For example, a higher acclimation temperature (57 days, 20 °C) increased the lytic activity of the total complement system in the plasma of rainbow trout [106]. Interestingly, the gene c3 (Complement C3), an important component of the alternative complement pathway was not significantly affected by exposure to elevated temperatures in the qPCR assessment. In accordance, we found that the same WN and WH fish did not show changes in hemolytic activity of the alternative pathway in the plasma 24 h after an intraperitoneal injection with bacterial and viral antigens in comparison to CT fish [47]. Similar findings of a partially activated classical pathway of the complement system without activating the alternative pathway, and without the induction of final cytotoxic activity, were reported in the liver of rainbow trout reared at 22 °C vs 14 °C [19]. Therefore, the alternative pathway, unlike the classical pathway, does not appear to have been activated under our experimental conditions. The classical pathway links innate and adaptive immunity as C1q binds to IgM molecules alone or aggregated on immune complexes [107]. Indeed, we found enrichment of ‘MHC complex binding’ processes and a trend for higher expression of mhcii (alias hla-a, Major Histocompatibility Complex II) in WH fish (WN: 1.0; WH: 1.5-fold, not qPCR validated). Thus, the applied thermal challenge may have activated adaptive immune responses, on which fish rely more at higher temperatures [108]. The GO/pathway term network also indicated that several HSP transcripts were associated with antigen processing and presentation (i.e., hla-a, hsp90aa1, hsp90ab1, hspa8, and tapbp) and the MHC class II binding complex (i.e., hsp90aa1 and hspa8). HSPs have been shown to act as chaperones for cytosolic peptides involved in MHC-driven antigen presentation to T-lymphocytes [84]. For instance, the cytosolic chaperone HSP90 is associated with peptide binding to MHC-class I and MHC-class II [109, 110]. Consequently, the observed up-regulation of hsp90aa1 transcripts in WN and WH fish as compared to CT fish, and its relation to ‘antigen processing pathways’ may have enhanced MHC peptide complex assembly for antibody production [84, 85, 109, 110].

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 (Matrix Metalloproteinase 9) (WN: 2.0; WH: 1.5-fold, not qPCR validated), hinting to potential tissue remodeling processes [111]. Matrix metalloproteinases (MMPs) are endopeptidases that cleave all structural elements of the ECM and are responsible for physiological and pathophysiological tissue remodelling [111, 112]. As such, wound healing, and pathological remodeling processes of damaged 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 a potential higher or more virulent pathogen abundance in a warmer aquatic environment [113, 114] as a ‘pre-adaptation’, or to initiate immune defence responses against invading pathogens or existing infections. Further, our data suggest that the combined stressors of high temperature and moderate hypoxia had a greater impact on the hepatic immunity of salmon. However, no significant clinical signs of infection or mortalities were recorded in this experiment [42], and when salmon of both warming scenarios were held at 20 °C for 4 weeks (WN and WH group) 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 [47]. Interestingly, WN and WH fish had a faster induction of the innate immune response in comparison to CT fish [47], suggesting that post-smolt Atlantic salmon exposed to 20 °C, under normoxia or hypoxia, were still immune competent. 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 [115, 116], were significantly increased (to ~ 30–40 ng mL^-1) in Atlantic salmon exposed to an incremental temperature elevation to 22 °C [117], with approx. 15% mortality at this temperature and 33% when temperatures reach 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. Nonetheless, additional research using live pathogen exposures in 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 results in accelerated mitochondrial ROS formation that can lead to oxidative stress and cellular damage [98, 118, 119]. Cellular oxidative stress due to prolonged hyperthermia can cause impaired mitochondria bioenergetics and structural alterations of cells and tissue [101, 120]. Tolerance to cellular oxidative stress is provided by an effective antioxidant system [119, 121], as well as an HSP response which attempts to maintain protein folding and mitochondrial integrity, and support cell function and survival [86, 87]. In our study, WN fish appeared to have a stimulated oxidative stress response (Theme ^#3) (i.e., ‘regulation of response to stress’ and ‘response to oxidative stress’) while WH fish showed up-regulated redox-related pathways (i.e., ‘antioxidant activity’, ‘positive regulation of release of cytochrome c from mitochondria’ and ‘thioredoxin reduction’). Hence, the activation of these particular oxidative stress responses (Theme ^#3), in addition to activated HSP, UPR, and ER-stress responses (Theme ^#1), 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 at 12 °C. 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, 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, 122, 123]. The expression of hif1α is considered to be a reliable hypoxic biomarker due to its up-regulation after acute hypoxia (i.e., hours) in several fish species such as Eurasian perch (Perca fluviatilis) [124], Atlantic croaker (Micropogonias undulatus) [125], and zebrafish (Danio rerio) [28]. However, results are not as consistent with regard to the effects of prolonged hypoxia. For example, while chronic hypoxia induced the up-regulation of hif1α transcripts in the ovaries of Atlantic croaker (21 days at 55% DO) [125] and in the liver of sea bass (Dicentrarchus labrax) (15 days at 51% DO) [126], it was down-regulated in the muscle and not affected in the liver of perch (15 days at 30% DO) [124]. Further, while it was unchanged in the liver of Atlantic salmon exposed to 4–5 mg O₂ L^-1 for 120 days at 12 °C, severe 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 et al. [119] 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). These author’s hypothesized that a more oxidized redox state during extreme heat could interfere (i.e., ‘switch-off’) with the hypoxic signaling response, and thus, prevent the complex HIF-1 induced physiological response. This hypothesis may be supported by the herein observed significant down-regulation of egln2 (alias phd1, Egl-9 Family Hypoxia Inducible Factor 2) in both warmed groups (qPCR validated) since it encodes for cellular oxygen sensor enzymes responsible for the post-translational regulation of HIF-1α proteins [127–129]. Three EGL-Nine homologs (EGLN1-3) regulate the abundance of HIF-1α proteins through proline hydroxylation and consequent proteasomal degradation [129], and were shown to be involved in signaling responses in the brain of large yellow croaker (Larimichthys crocea) exposed to acute hypoxia [130]. The effects of temperature stress on egl-9 homolog transcript expression have not been previously reported. Interestingly, we found a down-regulation of egln2 transcripts in both the WN and WH groups, and this may imply a temperature-dependent post-translational regulation of HIF1α in the liver of Atlantic salmon. In addition, the down-regulation of calm (alias cam, Calmodulin) in WN and WH fish as compared to CT fish (qPCR validated) suggests that there might be Ca²⁺/Calmodulin Kinase-dependent transcriptional regulation of HIF-1 [131]. The Ca²⁺/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. Taken together, 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. Clearly, 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 down-regulated in fish from the WN or 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 [132]. In line with our results, elevated temperatures (17–18 °C) caused 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₂ saturation) [133] and Atlantic croaker (four weeks at 1.7 mg DO L^-1) [134]. Further, the down-regulated gene ucp2 in WN and WH fish (qPCR validated), codes for Mitochondrial Uncoupling Protein 2 which uncouples oxidative phosphorylation from ATP synthesis resulting in energy dissipation [135]. Decreased expression of ucp2 transcript levels with increasing temperatures (15–25 °C) was previously reported in the gill and liver of pikeperch (Sander lucioperca) [92], and demonstrates that UCP2 has a thermogenic function [135]. For instance, 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 [136]. Under stress conditions, a 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 [137]. 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 [138]. Thus, alterations in mitochondrial function at high temperatures may be an important mechanism for thermal acclimation and thermal tolerance in this species, and further research is needed to unravel this mechanism.

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 signaling [139]. The up-regulation of prdx6 upon heat stress to protect the cell from oxidative stress has been reported in other marine animals [140, 141]. However, while Antarctic emerald rockcod (Trematomus bernacchii) exposed to warming temperatures had slightly increased expression of the prdx6-b paralog in the liver, the expression of the prdx6-a paralog was downregulated [141]. This suggests that prdx6 paralogs may have different temperature sensitivity characteristics, and that the transcriptional responses of prdx6 and its paralogs upon warming and hypoxia deserve further investigation.

Finally, the gene cirbp (Cold-Inducible RNA-Binding Protein) was equally down-regulated 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 [89, 142]. The cold-shock protein CIRBP acts as mRNA chaperone and 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 [142]. The lower expression of cirbp in this study is in accordance with several heat-stress studies on salmonid fishes [89].

Taken together, the up-regulation and enrichment of pathways related to oxidative stress (Theme ^#3), HSP-response (Theme ^#1) and apoptosis (Theme ^#4) in WN and WH fish suggest that the induction of antioxidant enzymes and redox pathways was an important line of defence against oxidative stress in these fish. Still, the oxidative stress response may have also been partly compromised at 20 °C because RNA and protein damage and/or degradation may not have been completely prevented. This would partly explain the down-regulation of several of our target genes related to oxidoreductase activity in the liver of post-smolt salmon, and suggests a decreased effectiveness of the antioxidant system under long-term stress conditions.

Cellular metabolism

The abiotic factors temperature and water oxygen level have a profound influence on the allocation of energy to maintenance versus growth in fish [3, 38, 143]. A reduction in metabolic processes can conserve energy during stressful conditions as imposed by thermal challenges [14, 17, 20, 22, 32, 99, 144, 145] or hypoxia [25, 28, 146, 147]. In this study, salmon exposed to the WN and WH conditions showed a suppression of genes related to a broad variety of cellular metabolic processes in the liver that were highly interconnected. 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’). For instance, the down-regulation 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 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 [145]. In the current study, the gene pdk3 was up-regulated in the liver of WN fish (qPCR validated) and this gene codes for pyruvate dehydrogenase kinase 3 (PDK3), which acts together with PDK1, PDK2 and PDK4 isoenzymes, to regulate glycolysis and glucose homeostasis under starvation [148].

In our study, salmon subjected to an incremental temperature increase to 20 °C and moderate hypoxia (WH group) 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 strong 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, 143]. For example, Atlantic salmon reared at 17–19 °C for 45 days showed reduced lipids and lower triacylglycerol (TAG) stores in the liver than in fish maintained at 13 °C, suggesting the reallocation and/ or depletion of endogenous lipid stores 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 [143]. 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, long-term exposure to high temperature and hypoxia may result in reduced lipid and fatty acid biosynthesis. In addition, rainbow trout exposed to an incremental temperature increase to 24 °C showed a down-regulation of hepatic genes related to energy metabolism in a temperature-dependent manner [20]. Thus, the down-regulation of pathways related to the aerobic metabolism of carbohydrates, proteins and fatty acids in the liver of stressed salmon may reflect a 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 hypoxic conditions, 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, as gene expression changes were highly associated with the growth performance (i.e., weight, length, and SGR) of the WH fish, and in this group we found positive correlations between growth parameters and the expression of 12 down-regulated genes (i.e. cirbp1, calm, egln2, hif1α, ucp2, gstt1, prdx6, rraga). Consequently, these changes in transcriptional molecular processes in WH fish could be connected to an impairment of their growth performance, and underline 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 which may have contributed to a metabolic suppression in our Atlantic 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 (qPCR validated), and this gene codes for DNA (cytosine-5)-methyltransferase, an enzyme essential for maintaining DNA methylation marks after mitosis [149]. 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., ~ 2,894 DEPs). Indeed, Beemelmanns et al. [150] found that the same treatments (i.e., WH and WH) affected the methylation of CpG sites of the herein 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 [150].

Biomarker genes

Transcriptomic techniques allow for the high-throughput 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 [89, 92, 151]. 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 [89, 92, 152]. The higher expression of serpinh1 and hsp90aa1 in different tissues during thermal stress is well documented in a number of different fish species [88, 89, 92, 152], and thus, our results further support their application as biomarkers for this stressor. HSP70 is usually considered to be a hallmark of the heat shock response in fishes [88, 89]. However, we observed 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 [153]. Ccomplementary to prior studies, we report marked dysregulation in the expression levels of 19 microarray-identified genes upon high temperature and hypoxia exposure that, in addition, showed associations with phenotypic characteristics. Hence, these genes may be useful not only as molecular biomarkers of thermal stress but as candidate genes for the development of thermal phenotype-relevant genomic markers [e.g., single nucleotide polymorphisms (SNPs) for marker-assisted selection of heat stress resistant broodstock], protein-based diagnostic assays (e.g., ELISA test) and for the detection of epigenetic markers (‘epimarkers’) that can predict thermotolerance [150].

To summarize, we have identified numerous transcriptional changes (~ 2,894 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 (~ 70% of air sat.), the latter realistically simulating summer conditions in salmon aquaculture sea cages [33, 34, 41]. These two treatments induced comparable biological processes with regards to the maintenance of cellular homeostasis (i.e., HSP, UPR, ER-stress response) and the apoptosis of damaged cells, and also stimulated immune gene expression responses (both innate and adaptive) to prepare for an anticipated higher, and possibly more virulent, pathogen load at warmer temperatures. In contrast, a large variety of general metabolic and proteolytic processes were reduced, and this likely reflects an extensive shift in liver metabolism that may have contributed to a metabolic suppression. This response may have been required for the reallocation of energy resources and to balance the energetically costly processes induced to maintain cellular homeostasis. Further, salmon exposed to either WN or WH conditions showed some signs of a compromised oxidative stress response at 20 °C, and this suggests a decreased effectiveness of the redox- and antioxidant system with long-term heat stress. Interestingly, the salmon that experienced the additional challenge of hypoxia had enhanced hepatic transcript expression related to the heat shock response and immunity as compared to WN fish, and the expression levels of 19 genes in the WH group were highly associated with reduced growth performance (i.e., weight, SGR, CF). These findings suggest that the combination of temperature stress and low oxygen levels had a greater impact on energy-demanding processes related to the maintenance of cell homeostasis (i.e., HSP, UPR, ER-stress immune response), and that this was balanced by a re-prioritization of energy expenditure in the liver that may have impaired the salmon’s growth and physiological performance. Consequently, these changes in transcriptional molecular processes in WH fish could be connected to phenotypic changes and underline the biological relevance of these hepatic transcriptional responses.

We expect that these results will be extremely valuable for the evaluation of thermal stress in fish in aquaculture and those exposed to climate-related challenges in the wild. Further, we anticipate that the 19 highly responsive and validated biomarker genes will contribute to the development of novel diagnostic tools such as protein-based assays (e.g., ELISA test), genomic markers (e.g., SNPs) and epigenetic markers (e.g., DNA methylation) to quantify and predict, more precisely, the stress levels and thermotolerance of wild and farmed fishes.

air sat: Air Saturation

CF: Condition Factor

CT: Control

C_T: Threshold Cycle

CP-1: First Component

CP-2: Second Component

CV: Coefficient of Variation

DEG: Differentially Expressed Gene

DEP: Differentially Expressed Probe

ECM: Extracellular Matrix

ER: Endoplasmic Reticulum

EST: Expressed Sequence Tag

FC: Fold Change

FDR: False Discovery Rate

GO: Gene Ontology

GOI: Gene of Interest

HSI: Hepato-Somatic Index

HSP: Heat Shock Protein

KEGG: Kyoto Encyclopedia of Genes and Genomes

lme: Linear Mixed Effect Model

NTC: Non-Template Control

PCA: Principal Component Analysis

PC: Principal Component

PC-1: First Principal Component

PC-2: Second Principal Component

qPCR: Quantitative PCR

ROS: Reactive Oxygen Species

RQ: Relative Quantity

RT: Reverse Transcription

RVM: Relative Ventricular Mass

SAM: Significance Analysis of Microarrays

SD: Standard Deviation

SEM: Standard Error of Mean

SGR: Specific Growth Rate

SSI: Spleen-Somatic Index

TAG: Triacylglycerol

UPR: Unfolded Protein Response

WH: Warm & Hypoxic

WN: Warm & Normoxic

Ethics approval and consent to participate

All experimental procedures described herein were approved by the Institutional Animal Care Committee of Memorial University (Protocol ^#16-90-KG) and followed guidelines set by the Canadian Council on Animal Care. This study also adheres to the ARRIVE Guidelines for reporting animal research [45], and a completed ARRIVE guidelines checklist is included as Additional file 1.

Consent for publication

This manuscript has not been published, nor is it under consideration elsewhere. All authors have agreed to the submission of this manuscript.

Availability of data and materials

The whole microarray dataset consisted of 17,072 detected probes (normalized log₂ ratios) and is accessible on-line on the NCBI’s Gene Expression Omnibus database through the GEO Series accession number GSE146471 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146471). The obtained threshold cycle (CT) values for samples from the current study (Fish-ID ^#73 to ^#96, sampling time point at 20°C/3-days) are accessible on-line at PANGAEA (https://doi.org/10.1594/PANGAEA.913696).

Competing interests

The authors declare that they have no known competing financial interests or personal relationships that influenced the work reported in this paper.

Funding

This study was conducted within the project ‘Mitigating the Impact of Climate-Related Challenges on Salmon Aquaculture (MICCSA)’ and was funded by the Atlantic Canada Opportunities Agency (781-9658-205222), Innovate NL (5404-1209-104) and Innovative PEI. QPCR primer and assay development was partly supported by the ‘Biomarker Platform for Commercial Aquaculture Feed Development’ project. This is a Genomic Applications Partnership Program (GAPP #6604) grant awarded to MLR and funded by the Government of Canada through Genome Canada and Genome Atlantic, and Cargill Innovation.

Authors' contributions

All authors contributed intellectually to this study. AB carried out the experiment and samplings, performed the microarray analysis, qPCR assay development, experimental qPCR measurements, statistical data analyses, data interpretation, and the writing of this manuscript. FZ contributed to the design and carried out the experiment and samplings. XX helped with the microarray design, analyses, and assisted with data analyses and interpretation. RS assisted with the experiment, sampling, and RNA sample preparation. MLR was involved in the experimental design, microarray experiment design, data analyses, and the interpretation of the results, and provided reagents. KG conceived the study, designed the experiment, provided reagents, and was involved with the interpretation of the results. All authors read, revised, edited, and approved the final manuscript.

Acknowledgments

The authors thank the staff of the Dr. Joe Brown Aquatic Research Building (JBARB, Memorial University, Canada) for assistance with fish husbandry. We are grateful to MacGregor Parent, Olufemi Ajiboye, Tasha Harrold and Gord Nash for assistance with sampling, and to Tasha Harrold for her help in managing the ‘Mitigating the Impact of Climate-Related Challenges on Salmon Aquaculture (MICCSA)’ project. We are thankful for the support of Jennifer Hall during primer quality testing and thank her and Albert Caballero-Solares for the provision of additional primers that were originally developed by them within a Genomic Applications Partnership Program (GAPP) project. We thank Olivia Roth and Thorsten Reusch for their collaboration in this project, and the use of their Fluidigm Biomark qPCR machine, and lab facilities and supplies at the Geomar (Kiel, Germany). We also thank Diana Gill for her help with sample shipment and technical assistance. Finally, we are grateful to Drs. Albert Caballero-Solares, Navaneethaiyer Umasuthan, and Khalil Eslamloo for many fruitful discussions about transcriptomic data analysis and interpretation, and to Drs. Surendra Kumar and Navaneethaiyer Umasuthan for providing an updated annotation of the 44K salmon microarray. We would also like to thank ACENET (https://www.ace-net.ca/) and Compute Canada (https://www.computecanada.ca/) for providing the computational resources to update the 44K microarray.

Wade NM, Clark TD, Maynard BT, Atherton S, Wilkinson RJ, Smullen RP, et al. Effects of an unprecedented summer heatwave on the growth performance, flesh colour and plasma biochemistry of marine cage-farmed Atlantic salmon (Salmo salar). J Therm Biol. 2019;80:64–74.
Motyka R, Norin T, Petersen LH, Huggett DB, Gamperl AK. Long-term hypoxia exposure alters the cardiorespiratory physiology of steelhead trout (Oncorhynchus mykiss), but does not affect their upper thermal tolerance. J Therm Biol. 2017;68:149–61. March 2016.
Vikeså V, Nankervis L, Hevrøy EM. Appetite, metabolism and growth regulation in Atlantic salmon (Salmo salar L.) exposed to hypoxia at elevated seawater temperature. Aquac Res. 2017;48:4086–101.
McBryan TL, Anttila K, Healy TM, Schulte PM. Responses to temperature and hypoxia as interacting stressors in fish: Implications for adaptation to environmental change. Integr Comp Biol. 2013;53:648–59.
Currie S, Schulte PM. Thermal stress. 4th editio. Boca Raton, Florida: CRC Press; 2014.
Kvamme BO, Gadan K, Finne-Fridell F, Niklasson L, Sundh H, Sundell K, et al. Modulation of innate immune responses in Atlantic salmon by chronic hypoxia-induced stress. Fish Shellfish Immunol. 2013;34:55–65.
Frölicher TL, Fischer EM, Gruber N. Marine heatwaves under global warming. Nature. 2018;560:360–4.
Oliver ECJ, Donat MG, Burrows MT, Moore PJ, Smale DA, Alexander LV, et al. Longer and more frequent marine heatwaves over the past century. Nat Commun. 2018;9:1324.
IPCC. 2019. Summary for Policymakers. In: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate [Pörtner H-O, Roberts DC, Masson-Delmotte V, Zhai P, Tignor M, Poloczanska E, Mintenbeck K, Alegría A, Nicolai M, Okem A, Petzold J, Rama B, Weyer NM, editors]. In press.
Breitburg D, Levin LA, Oschlies A, Grégoire M, Chavez FP, Conley DJ, et al. Declining oxygen in the global ocean and coastal waters. Science. 2018;359:80-. eaam7240.
Claret M, Galbraith ED, Palter JB, Bianchi D, Fennel K, Gilbert D, et al. Rapid coastal deoxygenation due to ocean circulation shift in the northwest Atlantic. Nat Clim Chang. 2018;8:868–72.
Kültz D. Molecular and evolutionary basis of the cellular stress response. Annu Rev Physiol. 2005;67:225–57.
Jeffries KM, Hinch SG, Sierocinski T, Clark TD, Eliason EJ, Donaldson MR, et al. Consequences of high temperatures and premature mortality on the transcriptome and blood physiology of wild adult sockeye salmon (Oncorhynchus nerka). Ecol Evol. 2012;2:1747–64.
Olsvik PA, Vikesa V, Lie KK, Hevroy EM. Transcriptional responses to temperature and low oxygen stress in Atlantic salmon studied with next-generation sequencing technology. BMC Genom. 2013;14:817.
Rebl A, Verleih M, Köbis JM, Kühn C, Wimmers K, Köllner B, et al. Transcriptome profiling of gill tissue in regionally bred and globally farmed rainbow trout strains reveals different strategies for coping with thermal stress. Mar Biotechnol. 2013;15:445–60.
Shi K-P, Dong S-L, Zhou Y-G, Li Y, Gao Q-F, Sun D-J. RNA-seq reveals temporal differences in the transcriptome response to acute heat stress in the Atlantic salmon (Salmo salar). Comp Biochem Physiol Part D Genomics Proteomics. 2019;30:169–78.
Jeffries KM, Hinch SG, Sierocinski T, Pavlidis P, Miller KM. Transcriptomic responses to high water temperature in two species of Pacific salmon. Evol Appl. 2014;7:286–300.
Tomalty KMH, Meek MH, Stephens MR, Rincón G, Fangue NA, May BP, et al. Transcriptional response to acute thermal exposure in juvenile chinook salmon determined by RNAsEq. Genes Genomes Genet. 2015;5:1335–49.
Oku H, Tokuda M, Matsunari H, Furuita H, Murashita K, Yamamoto T. Characterization of differentially expressed genes in liver in response to the rearing temperature of rainbow trout Oncorhynchus mykiss and their heritable differences. Fish Physiol Biochem. 2014;40:1757–69.
Li Y, Huang J, Liu Z, Zhou Y, Xia B, Wang Y, et al. Transcriptome analysis provides insights into hepatic responses to moderate heat stress in the rainbow trout (Oncorhynchus mykiss). Gene. 2017;619:1–9.
Liu S, Wang X, Sun F, Zhang J, Feng J, Liu H, et al. RNA-Seq reveals expression signatures of genes involved in oxygen transport, protein synthesis, folding, and degradation in response to heat stress in catfish. Physiol Genomics. 2013;45:462–76.
Jesus TF, Grosso AR, Almeida-Val VMF, Coelho MM. Transcriptome profiling of two Iberian freshwater fish exposed to thermal stress. J Therm Biol. 2016;55:54–61.
Guo L, Wang Y, Liang S, Lin G, Chen S, Yang G. Tissue-overlapping response of half-smooth tongue sole (Cynoglossus semilaevis) to thermostressing based on transcriptome profiles. Gene. 2016;586:97–104.
Richards JG. Chapter 10 Metabolic and Molecular Responses of Fish to Hypoxia. In: Fish physiology. Elsevier; 2009. p. 443–85.
Gracey AY, Troll JV, Somero GN. Hypoxia-induced gene expression profiling in the euryoxic fish Gillichthys mirabilis. Proc Natl Acad Sci. 2001;98:1993–8.
van der Meer DLM, van den Thillart GEEJM, Witte F, de Bakker MAG, Besser J, Richardson MK, et al. Gene expression profiling of the long-term adaptive response to hypoxia in the gills of adult zebrafish. Am J Physiol Integr Comp Physiol. 2005;289:R1512–9.
Qi D, Chao Y, Wu R, Xia M, Chen Q, Zheng Z. Transcriptome analysis provides insights into the adaptive responses to hypoxia of a schizothoracine fish (Gymnocypris eckloni). Front Physiol. 2018;9:1326.
Ton C, Stamatiou D, Liew C-C. Gene expression profile of zebrafish exposed to hypoxia during development. Physiol Genomics. 2003;13:97–106.
Wang Q-F, Shen W-L, Hou C-C, Liu C, Wu X-F, Zhu J-Q. Physiological responses and changes in gene expression in the large yellow croaker Larimichthys crocea following exposure to hypoxia. Chemosphere. 2017;169:418–27.
Chen B-X, Yi S-K, Wang W-F, He Y, Huang Y, Gao Z-X, et al. Transcriptome comparison reveals insights into muscle response to hypoxia in blunt snout bream (Megalobrama amblycephala). Gene. 2017;624:6–13.
Gong D, Xu L, Li W, Shang R, Chen J, Hu F, et al. Comparative analysis of liver transcriptomes associated with hypoxia tolerance in the gynogenetic blunt snout bream. Aquaculture. 2020;523:735163.
Prado-Lima M, Val AL. Transcriptomic characterization of tambaqui (Colossoma macropomum, Cuvier, 1818) exposed to three climate change scenarios. PLoS One. 2016;11:e0152366.
A Review of the 2019 Newfoundland and Labrador South Coast Cultured Atlantic Salmon Mortality Event St. John's

Burke H, Gardner I, Farrell AP. A Review of the 2019 Newfoundland and Labrador South Coast Cultured Atlantic Salmon Mortality Event St. John's, NL, Canada: Department of Fisheris and Land Resources, Government of Newfoundland and Labrador, Special Studies and Reports. 2020.
Burt K, Hamoutene D, Mabrouk G, Lang C, Puestow T, Drover D, et al. Environmental conditions and occurrence of hypoxia within production cages of Atlantic salmon on the south coast of Newfoundland. Aquac Res. 2012;43:607–20.
Gunderson AR, Armstrong EJ, Stillman JH. Multiple stressors in a changing world: The need for an improved perspective on physiological responses to the dynamic marine environment. Ann Rev Mar Sci. 2016;8:357–78.
Stehfest KM, Carter CG, McAllister JD, Ross JD, Semmens JM. Response of Atlantic salmon Salmo salar to temperature and dissolved oxygen extremes established using animal-borne environmental sensors. Sci Rep. 2017;7:4545.
Handeland SO, Imsland AK, Stefansson SO. The effect of temperature and fish size on growth, feed intake, food conversion efficiency and stomach evacuation rate of Atlantic salmon post-smolts. Aquaculture. 2008;283:36–42.
Hevrøy EM, Hunskår C, de Gelder S, Shimizu M, Waagbø R, Breck O, et al. GH-IGF system regulation of attenuated muscle growth and lipolysis in Atlantic salmon reared at elevated sea temperatures. J Comp Physiol B Biochem Syst Environ Physiol. 2013;183:243–59.
Solstorm D, Oldham T, Solstorm F, Klebert P, Stien LH, Vågseth T, et al. Dissolved oxygen variability in a commercial sea-cage exposes farmed Atlantic salmon to growth limiting conditions. Aquaculture. 2018;486:122–9.
Oldham T, Oppedal F, Dempster T. Cage size affects dissolved oxygen distribution in salmon aquaculture. Aquac Environ Interact. 2018;10:149–56.
Oppedal F, Vågseth T, Dempster T, Juell J-E, Johansson D. Fluctuating sea-cage environments modify the effects of stocking densities on production and welfare parameters of Atlantic salmon (Salmo salar L.). Aquaculture. 2011;315:361–8.
Gamperl AK, Ajiboye OO, Zanuzzo FS, Sandrelli RM, Peroni EFC, Beemelmanns A. The impacts of increasing temperature and moderate hypoxia on the production characteristics, cardiac morphology and haematology of Atlantic Salmon (Salmo salar). Aquaculture. 2020;519:734874.
Organic Chemical Toxicology of Fishes. Academic Press

Fish Physiology
Tierney KB, Farrell AP, Brauner CJ. Fish Physiology: Organic Chemical Toxicology of Fishes. Academic Press; 2013.
Jantzen SG, Sanderson DS, von Schalburg KR, Yasuike M, Marass F, Koop BF. A 44K microarray dataset of the changing transcriptome in developing Atlantic salmon (Salmo salar L.). BMC Res Notes. 2011;4:88.
Kilkenny C, Browne WJ, Cuthill IC, Emerson M, Altman DG. Improving bioscience research reporting: the ARRIVE guidelines for reporting animal research. PLoS Biol. 2010;8:e1000412.
Champely S, Ekstrom C, Dalgaard P, Gill J, Weibelzahl S, Anandkumar A, et al. Package ‘pwr.’ R Packag version. 2018;1.
Zanuzzo FS, Beemelmanns A, Hall JR, Rise ML, Gamperl AK. The innate immune response of Atlantic Salmon (Salmo salar) is not negatively affected by high temperature and moderate hypoxia. Front Immunol. 2020;11.
Xu Q, Feng CY, Hori TS, Plouffe DA, Buchanan JT, Rise ML. Family-specific differences in growth rate and hepatic gene expression in juvenile triploid growth hormone (GH) transgenic Atlantic salmon (Salmo salar). Comp Biochem Physiol Part D Genomics Proteomics. 2013;8:317–33.
Umasuthan N, Xue X, Caballero-Solares A, Kumar S, Westcott JD, Chen Z, et al. Transcriptomic Profiling in Fins of Atlantic Salmon Parasitized with Sea Lice: Evidence for an Early Imbalance Between Chalimus-Induced Immunomodulation and the Host’s Defense Response. Int J Mol Sci. 2020;21:2417.
Caballero-Solares A, Hall JR, Xue X, Eslamloo K, Taylor RG, Parrish CC, et al. The dietary replacement of marine ingredients by terrestrial animal and plant alternatives modulates the antiviral immune response of Atlantic salmon (Salmo salar). Fish Shellfish Immunol. 2017;64:24–38.
Xue X, Hixson SM, Hori TS, Booman M, Parrish CC, Anderson DM, et al. Atlantic salmon (Salmo salar) liver transcriptome response to diets containing Camelina sativa products. Comp Biochem Physiol Part D Genomics Proteomics. 2015;14:1–15.
Hori TS, Gamperl AK, Nash G, Booman M, Barat A, Rise ML. The impact of a moderate chronic temperature increase on spleen immune-relevant gene transcription depends on whether Atlantic cod (Gadus morhua) are stimulated with bacterial versus viral antigens. Genome. 2013;56:567–76.
Qiu W, Lee M-LT, Whitmore GA. Sample Size and Power Calculation in Microarray Studies Using the sizepower package. R Packag version. 2008;10:132.
Brazma A, Hingamp P, Quackenbush J, Sherlock G, Spellman P, Stoeckert C, et al. Minimum information about a microarray experiment (MIAME)-toward standards for microarray data. Nat Genet. 2001;29:365–71.
Yang YH, Paquet AC. Preprocessing two-color spotted arrays. In: Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Springer; 2005. p. 49–69.
Booman M, Xu Q, Rise ML. Evaluation of the impact of camelina oil-containing diets on the expression of genes involved in the innate anti-viral immune response in Atlantic cod (Gadus morhua). Fish Shellfish Immunol. 2014;41:52–63.
Bø TH, Dysvik B, Jonassen I. LSimpute: accurate estimation of missing values in microarray data with least squares methods. Nucleic Acids Res. 2004;32:e34–4.
Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci. 2001;98:5116–21.
Schwender H, Krause A, Ickstadt K. Identifying interesting genes with siggenes. Newsl R Proj. 2006;6:45.
Kolde R, Kolde MR. Package ‘pheatmap.’. R Packag. 2015;1:1–8.
Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25:1091–3.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.
Kanehisa M, Goto S, Kawashima S, Nakaya A. The KEGG databases at GenomeNet. Nucleic Acids Res. 2002;30:42–6.
Joshi-Tope G, Gillespie M, Vastrik I, D’Eustachio P, Schmidt E, de Bono B, et al. Reactome: A knowledgebase of biological pathways. Nucleic Acids Res. 2005;33 suppl_1:D428–32.
Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLoS One. 2011;6:e21800.
Caballero-Solares A, Xue X, Parrish CC, Foroutani MB, Taylor RG, Rise ML. Changes in the liver transcriptome of farmed Atlantic salmon (Salmo salar) fed experimental diets based on terrestrial alternatives to fish meal and fish oil. BMC Genom. 2018;19:796.
Soto-Dávila M, Valderrama K, Inkpen SM, Hall JR, Rise ML, Santander J. Effects of Vitamin D2 (Ergocalciferol) and D3 (Cholecalciferol) on Atlantic Salmon (Salmo salar) Primary Macrophage Immune Response to Aeromonas salmonicida subsp. salmonicida Infection. Front Immunol. 2020;10:3011.
Eslamloo K, Xue X, Hall JR, Smith NC, Caballero-Solares A, Parrish CC, et al. Transcriptome profiling of antiviral immune and dietary fatty acid dependent responses of Atlantic salmon macrophage-like cells. BMC Genom. 2017;18:706.
Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007;8:R19.
Beemelmanns A, Roth O. Biparental immune priming in the pipefish Syngnathus typhle. Zoology. 2016;119:262–72.
Bookout AL, Mangelsdorf DJ. Quantitative real-time PCR protocol for analysis of nuclear receptor signaling pathways. Nucl Recept Signal. 2003;1:nrs.01012.
R Core Team. R: A language and environment for statistical computing. 2018. http://www.r-project.org/.
Dray S, Dufour AB, Thioulouse J. ade4: Analysis of ecological data: Exploratory and Euclidean methods in environmental sciences. R package version 1.7-2. 2015.
Nguyen LH, Holmes S. Ten quick tips for effective dimensionality reduction. PLOS Comput Biol. 2019;15:e1006907.
Pinheiro J, Bates D, DebRoy S, Sarkar D. Linear and Nonlinear Mixed Effects Models. R Packag version. 2011;3:57.
Lenth RV. Least-squares means: the R package lsmeans. J Stat Sofw. 2016;69:1–33.
Beaton D, Chin Fatt CR, Abdi H. An ExPosition of multivariate analysis with the singular value decomposition in R. Comput Stat Data Anal. 2014;72:176–89.
Friendly M. Corrgrams: Exploratory displays for correlatigon matrices. Am Stat. 2002;56:316–24.
Reid G, Gurney-Smith H, Marcogliese D, Knowler D, Benfey T, Garber A, et al. Climate change and aquaculture: considering biological response and resources. Aquac Environ Interact. 2019;11:569–602.
Mohanty BP, Mahanty A, Mitra T, Parija SC, Mohanty S. Regulation of heat shock protein responses. In: Regulation of Heat Shock Protein Responses. Springer International Publishing; 2018. p. 71–94.
Hartl FU, Hayer-Hartl M. Protein folding. Molecular chaperones in the cytosol: From nascent chain to folded protein. Science. 2002;295:1852–8.
Parsell DA, Lindquist S. The function of heat-shock proteins in stress tolerance: degradation and reactivation of damaged proteins. Annu Rev Genet. 1993;27:437–96.
Srivastava P. Interaction of Heat Shock Proteins with Peptides and Antigen Presenting Cells: Chaperoning of the Innate and Adaptive Immune Responses. Annu Rev Immunol. 2002;20:395–425.
Callahan MK, Garg M, Srivastava PK. Heat-shock protein 90 associates with N-terminal extended peptides and is required for direct and indirect antigen presentation. Proc Natl Acad Sci. 2008;105:1662–7.
Lanneau D, Brunet M, Frisan E, Solary E, Fontenay M, Garrido C. Heat shock proteins: essential proteins for apoptosis regulation. J Cell Mol Med. 2008;12:743–61.
Ghosh S, Sarkar P, Basak P, Mahalanobish S, Sil PC. Role of heat shock proteins in oxidative stress and stress tolerance. In: Heat Shock Proteins and Stress. Cham: Springer International Publishing; 2018. pp. 109–26.
Roberts RJ, Agius C, Saliba C, Bossier P, Sung YY. Heat shock proteins (chaperones) in fish and shellfish and their potential role in relation to fish health: A review. J Fish Dis. 2010;33:789–801.
Akbarzadeh A, Günther OP, Houde AL, Li S, Ming TJ, Jeffries KM, et al. Developing specific molecular biomarkers for thermal stress in salmonids. BMC Genom. 2018;19:749.
Ishida Y, Nagata K. Hsp47 as a collagen-specific molecular chaperon. In: Methods in Enzymology. Academic Press; 2011. p. 167–82.
Wang Y, Liu Z, Li Z, Shi H, Kang Y, Wang J, et al. Effects of heat stress on respiratory burst, oxidative damage and SERPINH1 (HSP47) mRNA expression in rainbow trout Oncorhynchus mykiss. Fish Physiol Biochem. 2016;42:701–10.
Swirplies F, Wuertz S, Baßmann B, Orban A, Schäfer N, Brunner RM, et al. Identification of molecular stress indicators in pikeperch Sander lucioperca correlating with rising water temperatures. Aquaculture. 2019;501:260–71.
Jackson SE. Hsp90: Structure and Function. In: Molecular Chaperones. Springer Berlin Heidelberg; 2013. p. 155–240.
Buchner J, Jing L. Structure, function and regulation of the hsp90 machinery. Biomed J. 2013;36:106.
Kiang JG, Tsokos GC. Heat shock protein 70 kDa: molecular biology, biochemistry, and physiology. Pharmacol Ther. 1998;80:183–201.
Sitia R, Braakman I. Quality control in the endoplasmic reticulum protein factory. Nature. 2003;426:891–4.
D’Arcy MS. Cell death: a review of the major forms of apoptosis, necrosis and autophagy. Cell Biol Int. 2019;43:582–92.
Cheng CH, Yang FF, Liao SA, Miao YT, Ye CX, Wang AL, et al. High temperature induces apoptosis and oxidative stress in pufferfish (Takifugu obscurus) blood cells. J Therm Biol. 2015;53:172–9.
Li B, Sun S, Zhu J, Yanli S, Wuxiao Z, Ge X. Transcriptome profiling and histology changes in juvenile blunt snout bream (Megalobrama amblycephala) liver tissue in response to acute thermal stress. Genomics. 2019;111:242–50.
Poon WL, Hung CY, Nakano K, Randall DJ. An in vivo study of common carp (Cyprinus carpio L.) liver during prolonged hypoxia. Comp Biochem Physiol - Part D Genomics Proteomics. 2007;2:295–302.
Liu B, Xu P, Brown PB, Xie J, Ge X, Miao L, et al. The effect of hyperthermia on liver histology, oxidative stress and disease resistance of the Wuchang bream, Megalobrama amblycephala. Fish Shellfish Immunol. 2016;52:317–24.
Kruidering M, Evan G. Caspase-8 in Apoptosis: The Beginning of “The End”? IUBMB Life. 2000;50:85–90.
Weitzman JB, Fiette L, Matsuo K, Yaniv M. JunD protects cells from p53-dependent senescence and apoptosis. Mol Cell. 2000;6:1109–19.
Crespo-Sanjuán J, Zamora-Gonzalez N, Calvo-Nieves MD, Andres-Ledesma C. Apolipoprotein D. In: Advances in Lipoprotein Research. BoD–Books on Demand; 2017. p. 25.
Holland MCH, Lambris JD. The complement system in teleosts. Fish Shellfish Immunol. 2002;12:399–420.
Nikoskelainen S, Kjellsen O, Lilius EM, Schrøder MB. Respiratory burst activity of Atlantic cod (Gadus morhua L.) blood phagocytes differs markedly from that of rainbow trout. Fish Shellfish Immunol. 2006;21:199–208.
Nayak A, Pednekar L, Reid KB, Kishore U. Complement and non-complement activating functions of C1q: A prototypical innate immune molecule. Innate Immun. 2012;18:350–63.
Bowden TJ, Thompson KD, Morgan AL, Gratacap RML, Nikoskelainen S. Seasonal variation and the immune response: A fish perspective. Fish Shellfish Immunol. 2007;22:695–706.
Kunisawa J, Shastri N. Hsp90α chaperones large C-terminally extended proteolytic intermediates in the MHC class I antigen processing pathway. Immunity. 2006;24:523–34.
Rajagopal D, Bal V, Mayor S, George A, Rath S. A role for the Hsp90 molecular chaperone family in antigen presentation to T lymphocytes via major histocompatibility complex class II molecules. Eur J Immunol. 2006;36:828–41.
Yabluchanskiy A, Ma Y, Iyer RP, Hall ME, Lindsey ML. Matrix metalloproteinase-9: Many shades of function in cardiovascular disease. Physiology. 2013;28:391–403.
Caley MP, Martins VLC, O’Toole EA. Metalloproteinases and wound healing. Adv Wound Care. 2015;4:225–34.
Kimes NE, Grim CJ, Johnson WR, Hasan NA, Tall BD, Kothary MH, et al. Temperature regulation of virulence factors in the pathogen Vibrio coralliilyticus. ISME J. 2012;6:835–46.
Vezzulli L, Previati M, Pruzzo C, Marchese A, Bourne DG, Cerrano C. Vibrio infections triggering mass mortality events in a warming Mediterranean Sea. Environ Microbiol. 2010;12:2007–19.
Zanuzzo FS, Sabioni RE, Marzocchi-Machado CM, Urbinati EC. Modulation of stress and innate immune response by corticosteroids in pacu (Piaractus mesopotamicus). Comp Biochem Physiol Part A Mol Integr Physiol. 2019;231:39–48.
Pérez-Casanova JC, Rise ML, Dixon B, Afonso LOB, Hall JR, Johnson SC, et al. The immune and stress responses of Atlantic cod to long-term increases in water temperature. Fish Shellfish Immunol. 2008;24:600–9.
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 preparation.
Banh S, Wiens L, Sotiri E, Treberg JR. Mitochondrial reactive oxygen species production by fish muscle mitochondria: Potential role in acute heat-induced oxidative stress. Comp Biochem Physiol Part B Biochem Mol Biol. 2016;191:99–107.
Heise K. Oxidative stress during stressful heat exposure and recovery in the North Sea eelpout Zoarces viviparus L. J Exp Biol. 2006;209:353–63.
Belhadj Slimen I, Najar T, Ghram A, Dabbebi H, Ben Mrad M, Abdrabbah M. Reactive oxygen species, heat stress and oxidative-induced mitochondrial damage. A review. Int J Hyperth. 2014;30:513–23.
Machado C, Zaleski T, Rodrigues E, Carvalho CS, Cadena SMSC, Gozzi GJ, et al. Effect of temperature acclimation on the liver antioxidant defence system of the Antarctic nototheniids Notothenia coriiceps and Notothenia rossii. Comp Biochem Physiol Part - B Biochem Mol Biol. 2014;172–173:21–8.
Wenger RH. Cellular adaptation to hypoxia: O2-sensing protein hydroxylases, hypoxia‐inducible transcription factors, and O2 ‐regulated gene expression. FASEB J. 2002;16:1151–62.
Kilic M, Kasperczyk H, Fulda S, Debatin K-M. Role of hypoxia inducible factor-1 alpha in modulation of apoptosis resistance. Oncogene. 2007;26:2027–38.
Rimoldi S, Terova G, Ceccuzzi P, Marelli S, Antonini M, Saroglia M. HIF-1α mRNA levels in Eurasian perch (Perca fluviatilis) exposed to acute and chronic hypoxia. Mol Biol Rep. 2012;39:4009–15.
Rahman MS, Thomas P. Molecular and biochemical responses of hypoxia exposure in Atlantic croaker collected from hypoxic regions in the northern Gulf of Mexico. PLoS One. 2017;12:e0184341.
Terova G, Rimoldi S, Corà S, Bernardini G, Gornati R, Saroglia M. Acute and chronic hypoxia affects HIF-1α mRNA levels in sea bass (Dicentrarchus labrax). Aquaculture. 2008;279:150–9.
Semenza GL. HIF-1, O2, and the 3 PHDs. Cell. 2001;107:1–3.
Del Peso L, Castellanos MC, Temes E, Martín-Puig S, Cuevas Y, Olmos G, et al. The von Hippel Lindau/hypoxia-inducible factor (HIF) pathway regulates the transcription of the HIF-proline hydroxylase genes in response to low oxygen. J Biol Chem. 2003;278:48690–5.
Taylor MS. Characterization and comparative analysis of the EGLN gene family. Gene. 2001;275:125–32.
Liu W, Liu X, Wu C, Jiang L. Transcriptome analysis demonstrates that long noncoding RNA is involved in the hypoxic response in Larimichthys crocea. Fish Physiol Biochem. 2018;44:1333–47.
Yuan G, Nanduri J, Bhasker CR, Semenza GL, Prabhakar NR. Ca2+/calmodulin kinase-dependent activation of hypoxia inducible factor 1 transcriptional activity in cells subjected to intermittent hypoxia. J Biol Chem. 2005;280:4321–8.
Bistolas N, Wollenberger U, Jung C, Scheller FW. Cytochrome P450 biosensors—a review. Biosens Bioelectron. 2005;20:2408–23.
Olsvik PA. Effects of hypo- and hyperoxia on transcription levels of five stress genes and the glutathione system in liver of Atlantic cod Gadus morhua. J Exp Biol. 2006;209:2893–901.
Rahman MS, Thomas P. Effects of hypoxia exposure on hepatic cytochrome P450 1A (CYP1A) expression in Atlantic croaker: Molecular mechanisms of CYP1A down-regulation. PLoS One. 2012;7:e40825.
Brand MD, Esteves TC. Physiological functions of the mitochondrial uncoupling proteins UCP2 and UCP3. Cell Metab. 2005;2:85–93.
Bermejo-Nogales A, Calduch-Giner JA, Pérez-Sánchez J. Tissue-specific gene expression and functional regulation of uncoupling protein 2 (UCP2) by hypoxia and nutrient availability in gilthead sea bream (Sparus aurata): implications on the physiological significance of UCP1–3 variants. Fish Physiol Biochem. 2014;40:751–62.
Laskowski M, Augustynek B, Kulawiak B, Koprowski P, Bednarczyk P, Jarmuszkiewicz W, et al. What do we not know about mitochondrial potassium channels? Biochim Biophys Acta - Bioenerg. 2016;1857:1247–57.
Gerber L, Clow KA, Mark FC, Gamperl AK. Does Warm Acclimation Protect Atlantic Salmon (Salmo salar) Cardiac Mitochondrial Function Against the Negative Effects of Heat Stress? in preparation.
Arevalo J, Vázquez-Medina J. The role of Peroxiredoxin 6 in cell signaling. Antioxidants. 2018;7:172.
Park H, Ahn I-Y, Kim H, Cheon J, Kim M. Analysis of ESTs and expression of two peroxiredoxins in the thermally stressed Antarctic bivalve Laternula elliptica. Fish Shellfish Immunol. 2008;25:550–9.
Tolomeo AM, Carraro A, Bakiu R, Toppo S, Place SP, Ferro D, et al. Peroxiredoxin 6 from the Antarctic emerald rockcod: molecular characterization of its response to warming. J Comp Physiol B. 2016;186:59–71.
Zhong P, Huang H. Recent progress in the research of cold-inducible RNA-binding protein. Futur Sci OA. 2017;3:FSO246.
Kullgren A, Jutfelt F, Fontanillas R, Sundell K, Samuelsson L, Wiklander K, et al. The impact of temperature on the metabolome and endocrine metabolic signals in Atlantic salmon (Salmo salar). Comp Biochem Physiol Part A Mol Integr Physiol. 2013;164:44–53.
Chen X, Liu X, Li B, Zhang Q, Wang J, Zhang W, et al. Cold Inducible RNA Binding Protein Is Involved in Chronic Hypoxia Induced Neuron Apoptosis by Down-Regulating HIF-1α Expression and Regulated By microRNA-23a. Int J Biol Sci. 2017;13:518–31.
Fé-Gonçalves LM, Araújo JDA, dos Santos CHA, Val AL, Almeida-Val VMF. How will farmed populations of freshwater fish deal with the extreme climate scenario in 2100? Transcriptional responses of Colossoma macropomum from two Brazilian climate regions. J Therm Biol. 2020;89:102487.
Zhong X-P, Wang D, Zhang Y-B, Gui J-F. Identification and characterization of hypoxia-induced genes in Carassius auratus blastulae embryonic cells using suppression subtractive hybridization. Comp Biochem Physiol Part B Biochem Mol Biol. 2009;152:161–70.
Li M, Wang X, Qi C, Li E, Du Z, Qin JG, et al. Metabolic response of Nile tilapia (Oreochromis niloticus) to acute and chronic hypoxia stress. Aquaculture. 2018;495:187–95.
Kuntz MJ, Harris RA. Pyruvate Dehydrogenase Kinase. Encycl Signal Mol. 2018;:1–9.
Edwards JR, Yarychkivska O, Boulard M, Bestor TH. DNA methylation and DNA methyltransferases. Epigenetics Chromatin. 2017;10:23.
Beemelmanns A, Ribas L, Anastasiadi D, Moraleda-Prados J, Zanuzzo SF, Rise ML, Gamperl AK. DNA methylation dynamics in Atlantic salmon (Salmo salar) challenged with high temperature and moderate hypoxia. submitted. 2020.
Eissa N, Wang H-P. Transcriptional stress responses to environmental and husbandry stressors in aquaculture species. Rev Aquac. 2016;8:61–88.
Houde ALS, Akbarzadeh A, Günther OP, Li S, Patterson DA, Farrell AP, et al. Salmonid gene expression biomarkers indicative of physiological responses to changes in salinity and temperature, but not dissolved oxygen. J Exp Biol. 2019;222:jeb198036.
Webb D, Gagnon MM. The value of stress protein 70 as an environmental biomarker of fish health under field conditions. Environ Toxicol. 2009;24:287–95.

Download PDF

Journal Publication

published 12 Apr, 2021

Read the published version in BMC Genomics →

Editorial decision: Minor revision
16 Jan, 2021
Review #2 received at journal
13 Jan, 2021
Review #1 received at journal
26 Dec, 2020
Reviewer #2 agreed at journal
13 Dec, 2020
Reviewer #1 agreed at journal
17 Aug, 2020
Reviewers invited by journal
07 Aug, 2020
Editor assigned by journal
09 Jul, 2020
Submission checks completed at journal
08 Jul, 2020
Editor invited by journal
08 Jul, 2020

You are reading this latest preprint version

The Transcriptomic Responses of Atlantic Salmon (Salmo salar) to High Temperature Stress Alone, and in Combination with Moderate Hypoxia

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Methods

Animal husbandry / initial holding

Experimental protocol

RNA extraction, DNase treatment, and column purification

Microarray protocol

Microarray experimental design and hybridization

Microarray data acquisition

Microarray data analysis

Functionally organized GO/pathway term network analysis

qPCR protocol

Gene selection and primer design

cDNA synthesis and primer efficiencies

Selection of normalizer genes

qPCR measurements (Fluidigm Biomark™)

qPCR data acquisition

Statistical analyses

Principal Component Analysis (PCA)

Gene-by-gene analysis

Correlation analyses

Results

Significance Analysis of Microarrays (SAM)

GO/pathway term network analysis

GO/pathway term networks associated with up-regulated DEGs

Heat Shock Response (#1)

Cellular Stress (#2) and Oxidative Stress (#3)

Apoptosis (#4)

Immune Response (#5)

Protein Processing & Localization (#6) and Transcription (#7)

GO/pathway term networks associated with down-regulated DEGs

Proteolysis (#8), Catabolic Processes (#9), and Cellular Metabolic Processes (#10)

qPCR validation of microarray approach

Gene-by-gene analysis

Genes related to Heat Shock Response (#1)

Genes related to Apoptosis (#4)

Genes related to Immune Response (#5)

PCA analysis on 24 stress-related (#1–4) and 14 immune-related (#5) genes

Correlation between gene expression and phenotypic traits

Correlation analysis for WN fish

Correlation analysis for WH fish

Discussion

Heat shock response, unfolded protein response and endoplasmic reticulum stress

Apoptosis and programmed cell death

Immune response

Oxidative stress

Cellular metabolism

Transcriptional regulation and epigenetic mechanisms

Biomarker genes

Conclusions

List Of Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Competing interests

Funding

Authors' contributions

Acknowledgments

References

Supplementary Files

Status:

Journal Publication

Version 1

qPCR measurements (Fluidigm Biomark^™)

Heat Shock Response (^#1)

Cellular Stress (^#2) and Oxidative Stress (^#3)

Apoptosis (^#4)

Immune Response (^#5)

Protein Processing & Localization (^#6) and Transcription (^#7)

Proteolysis (^#8), Catabolic Processes (^#9), and Cellular Metabolic Processes (^#10)

Genes related to Heat Shock Response (^#1)

Genes related to Apoptosis (^#4)

Genes related to Immune Response (^#5)

PCA analysis on 24 stress-related (^#1–4) and 14 immune-related (^#5) genes