Low Oxygen Stress During Early Development Influences Regulation of Hypoxia-Response Genes in Farmed Atlantic Salmon (Salmo salar)

Survival and growth of developing salmonids are negatively affected by low oxygen levels within gravel nests in natural streams, and hypoxic stress is often experienced by farmed Atlantic salmon (Salmo salar) within hatcheries. Exposure to hypoxia during early development may have long-lasting effects by altering epigenetic marks and gene expression in oxygen regulatory pathways. Here, we examine the transcriptomic response to low dissolved oxygen (DO) in post-hatch salmon reared continuously in 30%, 60% or 100% DO from fertilization until start of feeding. RNA sequencing revealed multiple differentially expressed genes, including oxygen transporting hemoglobin embryonic α subunit (hbae) and EGLN3 family hypoxia-inducible factor 3 (egln3) which regulates the stability of hypoxia inducible factor 1α (HIF-1α). Both hbae and egln3 displayed expression levels inversely correlated to oxygen concentration, and DNA methylation patterns within the egln3 promoter were negatively associated with the transcript levels. These results suggest that epigenetic processes are influenced by low oxygen levels during early development in Atlantic salmon to upregulate hypoxia-response genes.

. When such changes increase the chance of survival, such a process is known as adaptive plasticity (Ghalambor, McKay, Carroll, & Reznick 2007). A growing focus in aquaculture research has been to limit the physiological effect of environmental stress, in order to improve fish health and production efficiency. This has led to an interest in the role of epigenetics in gene-environment interactions in aquaculture. Atlantic salmon exposed to environmental stress during early developmental stages show modified transcription and methylation associated with body growth and immune response (Moghadam et al. 2017;Uren-Webster et al. 2018). Parental hypoxic exposure in zebrafish (Danio rerio) increases hypoxia tolerance in offspring (Ho and Burggren 2012), rearing chinook salmon (O. tshawytscha) embryos in hypoxia improves fry tolerance to acute hypoxia stress (Del Rio, Davis, Fangue, & Todgham 2019) and the hatchery environment induces epigenetic modifications in coho salmon (O. kisutch) (Le Luyer et al., 2017) and steelhead trout (O. mykiss) (Gavery et al., 2019). Early exposure to chronic hypoxia in European sea bass (Dicentrarchus labrax) had short-and long-term effects on expression of hemoglobins, several hypoxia inducible factor 1 (HIF-1) related genes and egl-9 family hypoxia-inducible factor 3 (egln3), regulating the stability of the HIF-1a subunit (Vanderplancke et al. 2015;Cadiz et al. 2017).
In the current study, we investigated the epigenetic regulation of the hypoxic response in Atlantic salmon by analyzing the transcriptome and methylome in alevins after embryonic exposure to various chronic hypoxic treatments. This is part of a broader investigation of the potential to condition Atlantic salmon during early life stages for tolerance to hypoxia.

Experimental set up and sampling
Atlantic salmon eggs and milt were provided by AquaGen AS (Trondheim, Norway) and fertilized according to standard hatchery "dry fertilisation" protocols. Eggs from one female, including ovarian fluid, and milt from one male (2ml per liter egg) were gently mixed for 2 min prior to adding water. Milt and ovarian fluid were gently washed away, and eggs were left to swell and harden for 2 hr, then moved to incubators. Eggs that were unfertilized, underdeveloped, dead or ceased developing at the blastula stage, were visually identified and removed.
Immediately following fertilization, eggs were divided into three oxygen treatment groups i) 30% dissolved oxygen (30% O 2 ) ii) 60% dissolved oxygen (60% O 2 ) and iii) 100% dissolved oxygen (100% O 2 ) and incubated in circular upwelling incubators (21cm diameter) at 7.3°60.3°, at 7.2°60.4°and at 6.7°60.3, respectively. Dissolved oxygen concentrations of 30% and 60% were within the range that salmon may encounter in the sea-cage environment (26-90%; Stehfest et al. 2017;Solstorm et al. 2018). The 100% dissolved oxygen condition was used as a control. The 30% dissolved oxygen concentration was achieved by degasification using pure nitrogen gas when trickling water through an approximately 1m 3 box filled with plastic media to increase the surface area. The 60% saturation treatment was achieved by mixing untreated water and 30% saturated water. The process of degasification caused the slightly higher temperature for lower oxygen treatments as detailed above. Water flow was checked weekly and adjusted to approximately the same rate in all incubators (variance 6 0.1 L/min).
Each treatment group was kept in triplicate tanks with approximately 1300 eggs per replicate. Eggs rested upon perforated plates (4mm diameter holes) within the incubators. Once hatched, larvae swam through the openings and rested on the plastic substrate (Figure 1).
At several timepoints throughout the alevin stage ( Figure S1), individuals were randomly sampled and directly euthanised using an overdose of Benzocain (Benzoak Vet.). Whole alevin samples without yolk sac were placed in RNA-later (Ambion) and kept at 4°overnight before storage at -20°. A Kruskal-Wallis test was conducted in R (R Core Team 2019) to assess if weight differences between treatment groups was significant (n = 120, P , 0.05). Ethical guidelines supplied by the Norwegian Food Safety authority were followed. The guidelines state that experiments using fish prior to start-feeding are not subject to ethics approval.

RNA and DNA extraction
Total RNA and DNA was isolated and purified from whole alevins without yolk sac using AllPrep DNA/RNA/miRNA Universal Kit (Qiagen) according to the manufacturer's protocol. RNA quantity and quality were determined using nanodrop (Thermo Scientific) and bioanalyzer (Agilent) instruments. Isolated RNA with a 260/280 ratio .1.9 and RNA integrity number (RIN) .8, and DNA with a 260/280 ratio .1.8 was considered to be pure. Isolated DNA and RNA was stored until further analysis at -20°and -80°, respectively.
RNA sequencing, alignment and differential expression analysis At 925 degree days (DD), just before start of feeding, twelve alevins (4 per replicate) from each of the three oxygen treatment groups were used for RNA sequencing. We focused on this point of development because larvae up to this point had been solely dependent on their endogenous nutrient supply (yolk), and we were concerned that after the salmon lose the yolk sac gene expression patterns might be influenced by exogenous feeding. Due to the impact of low oxygen saturation on development, individuals reared in 30% O 2 were visibly less developed at 925DD than those at 60% O 2 and 100% O 2 , marked by less absorption of the yolk sac. To minimize developmental stage effects, an additional 12 individuals were randomly sampled from the 30% O 2 treatment group at 990DD, which matched the same developmental stage (end of yolk sac absorption) as the other groups.
Illumina sequencing of the 48 RNA samples was conducted by the Norwegian Sequencing Centre (Oslo, Norway). Libraries were prepared using a Strand-specific TruSeq RNA-seq library preparation kit (Illumina), and paired-ends sequenced using a HiSeq 3/4000 Genome Analyzer (Illumina), in 4 lanes with a read length of 75bp.
Leading and trailing low quality or N bases (quality score ,10) were removed from the sequence using Trimmomatic (Bolger et al. 2014). The reads were scanned using a 4-base window, and the window was cut when the average quality per base fell below 15. Reads less than 36 bases in length were removed from the analysis.
Reads were mapped to the Atlantic salmon genome assembly (GCA000233375.4 ICSASG_v2.) (Lien et al. 2016) using TopHat v.2.0.8b with default parameters other than b2 optionsensitive Langmead et al. 2009). SamTools (Li et al. 2009) was used to sort the alignments and HtSeq-count was used to count how many reads mapped to each interval on each chromosome, skipping reads with alignment quality below 10 (Anders et al. 2015).
The DESeq2 package in R was used to quantify differential gene expression between oxygen treatments using a negative binomial linear model with default parameters (Love et al. 2014). Genes from the DESeq2 analysis with adjusted p-values , 0.05 were considered as differentially expressed. A complete linkage method was used to perform a hierarchical cluster analysis using the heatmap.2 and hclust packages in R.

Enrichment analysis of Gene Ontology terms
The R package topGO (Alexa et al. 2006) was used to perform an enrichment analysis for Gene Ontology (GO) terms (Ashburner et al. 2000; The Gene Ontology Consortium 2019) using the list of differentially expressed genes detected between oxygen treatment groups against all genes in the annotated reference genome. Significance level was established using Fisher's exact test (P , 0.05). GO annotations for Salmo salar reference sequence ICSASG_v2 were downloaded from GitHub (https://rdrr.io/github/FabianGrammes/ Ssa.RefSeq.db/).

qPCR gene expression analysis
Quantitative real-time PCR (qPCR) was used to investigate relative expression of important genes identified by RNA-seq at two additional sampling timepoints: 1) the point of hatching which occurred at around 560DD (and is therefore denoted here as 560DD) and 2) 700DD.
Total RNA was extracted from twelve individuals sampled per treatment group and cDNA was synthesized using the High-Capacity RNA-to-cDNA Kit (Applied Biosystems) according to manufacturer's protocol, using 200ng of total RNA. Significantly differentially expressed genes with putative functions affecting hypoxic response were targeted for qPCR. Target genes were egln3, egln2, egln1_L and hemoglobin embryonic a (hbae). Egln1_L and egln2 were not among the differentially expressed genes according to RNAseq analysis but were included as qPCR targets for comparison. Primers were designed using Primer Express 3 software (Life Technologies) (Table  S1). To ensure the specificity of primers, paralogous gene sequences were aligned using MUSCLE (https://www.ebi.ac.uk/Tools/msa/muscle/), viewed using GenDoc Multiple Sequence Alignment Editor & Shading Utility Version 2.7 (https://genedoc.software.informer.com/ 2.7/), primers designed in areas of sequence difference between paralogs and primer sequences checked for secondary matches using BLAST at NCBI. Once the primers were ordered, checks were made for primer dimers, secondary amplification products and secondary melting curve peaks using standard methods.
A twofold standard curve of cDNA was prepared, combing 12 samples across all treatments and creating an eight-point dilution series. This was done for each primer set to determine amplification efficiency and a melting curve analysis was used to check for primerdimers and untargeted amplification products. No template controls, and a control to confirm the absence of genomic DNA, were run. The qPCR was run in duplicates using QuantStudio 5 (Thermo Fisher). Each well contained Power SYBR Green (Thermo Fisher), 300nM final concentration of each primer, 7mL of diluted cDNA and nuclease free water (Ambion) to a final reaction volume of 20 mL.
Quant Studio Design & Analysis Software (Thermo Fisher), was used to collect and analyze all data. For each gene, the sample with the Figure 1 Schematic of upwelling incubator setup showing water flow (blue arrows) and oxygen saturation. Eggs rested upon perforated plates (blue). Once hatched, larvae swam through the openings and lay on a plastic substrate (gray). After hatching the perforated plate was removed.

Figure 2
Weight (g) of Atlantic salmon alevins (without yolk sac) sampled for the various oxygen treatments and timepoints throughout the experiment. The 990DD stage of the 30% O 2 group was included in the analyses due to slowed development as explained in M&M. Outliers are plotted as individual points. Letters indicate significant weight differences between treatment groups at each point of time (P , 0.05, Kruskal-Wallis test). lowest expression (highest mean CT value) across all time points was used as what is denoted as a calibrator by (Pfaffl et al. 2004). Reference genes were used as internal controls to normalize expression data and remove non-biological variation. The reference genes ef1a, 18S and b-actin were verified using BestKeeper software (Pfaffl et al. 2004) and the geometric mean of reference gene CT values was used to calculate relative expression of target genes (Vandesompele et al. 2002). Significance was tested using a Kruskal-Wallis test.

Analysis of DNA methylation in egln3 promoter
Twelve samples from each oxygen treatment group at the 700DD sampling timepoint were used for analysis of DNA methylation percentage on a PyroMark Q24 advanced system (Qiagen) (n = 36). This sampling time point was chosen because we hypothesized that most of the effects of treatment on DNA methylation would have occurred by this stage and that levels of DNA methylation would persist until differential gene expression was assessed.
Bisulfite conversion of DNA was completed using a Epitect Fast Bisulfite Conversion Kit (Qiagen) according to the protocol. Two pyrosequencing assays were designed using a PyroMark Assay Design 2.0 (Qiagen) to cover five of the eight CpG sites in the Atlantic salmon egln3 putative promoter region (Table S2). Primers could not be adequately designed to cover CpG sites 1-3. PCR was carried out using the PyroMark PCR Kit (Qiagen) according to the protocol with Mg 2+ added to the reaction to improve the signal. Singular PCR products were confirmed on a 1% agarose gel. PyroMark Q24 Advanced CpG Reagents (Qiagen) were used along with Streptavidin Sepharose High Performance Beads (GE Healthcare). The percentage of methylation per CpG site was obtained by pyrosequencing performed on a PyroMark Q24 pyrosequencer following manufacturer's protocol. Briefly, samples were first washed on a PyroMark Q24 Vacuum Workstation (Qiagen) using 70% ethanol, denatured using PyroMark Denaturation Solution (Qiagen) and washed in PyroMark Washing Buffer (Qiagen). Samples were then annealed with the sequencing primer in a PyroMark Q24 Plate (Qiagen), which was thereafter warmed for 5 min at 80°using a pre-heated PyroMark Q24 plate holder. Within 30 s, the heated plate containing single stranded template DNA with sequencing primer was transferred into the PyroMark Q24 Instrument and the run was started to collect the methylation data. Finally, PyroMark Q24 Advanced software was used to visualize and collect the methylation percentage per site and assess the quality of the data following manufacturer's guidelines.
The Kruskal-Wallis test was used in R to determine if any of the five CpG sites in the egln3 promoter showed significantly different methylation between the treatment groups (n = 12 per treatment, P , 0.05). Association between DNA methylation changes in the egln3 promoter and relative abundance of egln3 transcripts was examined using DNA and RNA extracted from the same individuals.
Methylation and expression data for the 36 individuals was collected and the cor.test package in R was used to assess correlation between the data, based on Pearson's correlation coefficient, following a t-distribution with n-2 degrees of freedom and a significance threshold of 0.05.

Data availability
RNA sequence data for this study has been submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra under accession number PRJNA609825). File S1 contains results of differential gene expression analysis. File S2 contains all other supplementary figures and tables. Supplemental material available at figshare: https:// doi.org/10.25387/g3.12587207.

Alevin weight
The weight of alevins was positively associated with dissolved oxygen across all timepoints (560DD, 700DD and 925DD, Figure 2). The weight disparity between the 30% O 2 group and the other groups at 925DD (P , 0.05), in addition to the presence of a yolk sac in only the 30% O 2 group at 925DD, indicate that this group was relatively underdeveloped at this timepoint. This was supported by the Principal Component Analysis (PCA) of the transcriptomes separating the 925DD 30% O 2 cluster from the 925DD 60% O 2 , 925DD 100% and 990 30% O 2 clusters (Figure 3, principal component 1). To reduce any potentially confounding effects of developmental stage we have focused mainly on the comparison of 30% O 2 990DD, 60% O 2 925DD and 100% O 2 925DD treated fish.
Differential gene expression associated with oxygen treatment Samples clustered according to oxygen treatment in the PCA with 30% O 2 , 60% O 2 and 100% O 2 treated fish distributed from left to right along the x-axis of the PCA plot (Figure 3). In this way PC1, which explained most (22%) of the variance, was positively associated with treatment oxygen level. The factor influencing PC2 which explained 18% of the variance was unknown. Although the unknown influencing factor on PC2 might reduce the power to detect differential expression associated with treatment, PC2 is clearly not confounded with treatment, and therefore should not be causing  false positive results in the analysis of treatment effects. Analysis of RNAseq data (File S1) identified 39 genes that were systematically differentially expressed when comparing expression between the 900DD-30% O 2 and 925DD-60% O 2 treatments, 900DD-30% O 2 and 925DD-100% O 2 treatments and the 60% O 2 and 100% O 2 treatments (Figure 4). The expression of most of the 39 genes, including egln3, hbae, cytochrome c oxidase-like cox4 and cox8, and DNA damage-inducible transcript 4-like (ddit4), were upregulated with reduced % O 2 treatments in a dose-dependent manner ( Figure 5, Table S3). Comparison of the two most extreme treatments (30% O 2 vs. 100% O 2 ) identified the largest numbers of differentially expressed genes (7738) and many of these genes (3032) were also differentially expressed in the other treatment comparisons. A lower number of differentially expressed genes (174) were detected for the 30% O 2 to 60% O 2 treatment comparison than for the other comparisons.
Gene Ontology terms that were significantly over-represented when comparing 60% O 2 and 100% O 2 treatments included "structural constituent of ribosome", "oxidoreductase activity, acting on heme group of donors", "heme-copper terminal oxidase activity" and "cytochrome-c oxidase activity", including various cytochrome c oxidase subunit genes in the respiratory electron transport chain. Many of these terms contain genes typically involved in redox reactions (Table S4). The GO-terms were slightly different with comparison of 30% O 2 to 60% O 2 treatments, and additional terms such as "iron ion binding", "cofactor binding", "molecular carrier activity", "oxygen binding" and "heme binding" were over-represented with this comparison (Table S4). A mixture of the same terms described for the two comparisons above were over-represented when comparing 30% O 2 and 100% O 2 treatments. Many new terms not observed in the above three comparisons, such as "catalytic activity" (including 2-oxoisovalerate dehydrogenase subunits), "isomerase activity" (including disulphide isomerases and enolase), "cofactor binding" and "translation factor activity" were over-represented when comparing 990DD and 925DD developmental stages for the 30% O 2 treatment (Table S4). Some terms in common with the O 2 treatment comparisons above were also over-represented in the comparison of 30% O 2 treatment at 990DD and 925DD, such as "oxidoreductase activity", "coenzyme binding" and "vitamin binding".
These terms therefore represent genes that responded differently to the 30% O 2 treatment depending on developmental stage.
We decided to focus on the 39 genes that were differentially expressed between all the O 2 treatment comparisons, in particular the hypoxia-related genes egln3 and hbae ( Figure 5). Egln3 (XM_ 014194431.1) had the most significant difference (lowest p-value) in transcript abundance between the three groups (with the highest transcript abundance in the 30% O 2 group, followed by 60% O 2 , and then 100% O 2 groups (Figure 6a). Hbae transcript levels (NM_ 001140924.1) showed a similar trend to that of egln3 (Figure 6b). These trends also held at the earlier developmental stage of 30% O 2 925DD, but the egln3 transcript levels were slightly higher in 30% O 2 925DD individuals compared to 30% O 2 990DD individuals (nonsignificant, Figure S2).
The relative expression levels of egln3, egln2, egln1_L and hbae were quantified by qPCR in whole alevins at earlier stages 560DD (hatching) and 700DD. Egln3 expression was negatively associated Figure 6 Transcript abundance (normalized counts, scaled to account for sequencing depth and RNA composition according to the method of Anders and Huber 2010) for A. egln3 and B. hbae in Atlantic salmon reared in 30% O 2 , 60% O 2 or 100% O 2 at 925DD (990DD for 30% O 2 group). Outliers are plotted as individual points. Letters indicate significant differential gene expression between treatment groups (P , 0.05).
with the O 2 levels at 560 and 700DD (Figure 7a). Egln3 expression was significantly upregulated with the 30% DO treatment at 990DD compared to 60 and 100% treatments at 925DD (Figure 7a).
The relative expression levels of egln2 and egln1_L followed a similar pattern to egln3 at 560DD and 700DD (Figure 7b and c). hbae expression was reduced as development proceeded and showed significant upregulation with 30% DO treatment at 700DD and 990DD (Figure 7d).

DNA methylation in egln3 promoter
The impact of low oxygen treatments on the DNA methylation of the CpG sites numbered 4-8 in the egln3 promoter were investigated at 700DD. Methylation decreased at site 4 as oxygen treatment increased, while at site 5 the greatest average methylation was observed in the 60% O 2 treatment group, with the 30% O 2 and 100% O 2 groups having a lower, but similar average methylation. A difference in methylation between treatment groups was detected at sites 6 and 8 (P , 0.05, Table S3, Figure 8).
The correlation between egln3 expression and DNA methylation at promoter CpG site 6 (r = -0.497, P , 0.01) and CpG site 8 (r = -0.494, P , 0.01) was negative and moderate in strength (Figure 9). In addition, the correlation between expression and methylation at CpG sites 6 and 8 combined was significant (r = -0.56, P , 0.01) and greater than that at sites 6 and 8 alone.

DISCUSSION
Environmental stresses such as low dissolved oxygen during early life can influence phenotype and can potentially affect the animal's response to this stressor at later stages (Luo and Wang 2018;Figure 7 Relative gene expression (qPCR) of three Atlantic salmon egln paralogs and hbae in alevins at 560DD and 700DD. Egln3 (a), egln2 (b), egln1L (c) and hbae (d). Outliers are plotted as individual points. Letters indicate significant differential gene expression between treatment groups at each point of time (P , 0.05, Kruskal-Wallis test). Treatments for time points 925DD (60% and 100% O 2 ) and 990DD (30% O 2 ) were also compared.  Nanduri et al. 2017;Hancock et al. 2015). Our results show that low oxygen treatments during early development of farmed Atlantic salmon have a significant effect on transcript expression of hypoxicresponse genes. The negative correlation between egln3 expression and DNA methylation in two CpG sites of the egln3 promoter suggests that hypoxic treatments trigger changes in DNA methylation within the promoter to upregulate egln3. Accordingly, CpG methylation of the human egln3 homolog was found to prevent hypoxia inducible factor (HIF)-1a degradation in plasma cell neoplasia leading to the activation of HIF target genes (Shah et al. 2007;Hatzimichael et al. 2010). Acute low oxygen treatments (1% O 2 for 1-2 h) have been found to induce egln3 expression during the first 24 h of development in zebrafish (Manchenkov et al. 2015) at a stage when the HIF pathway has been considered inactive (Mendelsohn and Gitlin 2008;Robertson et al. 2014). Applying extreme acute hypoxic treatments (0.5% and 5% dissolved oxygen for 4 h) to zebrafish embryos induced the HIF-1 oxygen response pathway and improved hypoxia tolerance in larvae, however this improved tolerance was not observed to persist in juvenile and adult developmental stages (Robertson et al. 2014). Together these studies indicate that the upregulation of egln family genes induced during treatment of the embryo may persist after the alevin stage of development and suggest that CpG methylation of the egln3 promoter contributes to regulate the expression level of the gene. Methylation analysis of juvenile and adult salmon will help to assess the potential for longer-lasting gene expression and methylation changes and determine whether such treatments could condition the salmon epigenome for resilience to hypoxia.
The EGLN enzymes, or mammalian PHD homologs, regulate the stability of HIF during normoxic conditions by targeting HIF-1a for degradation (Appelhoff et al. 2004;Maxwell et al. 1999;Ivan et al. 2001;Jaakkola et al. 2001;Smith et al. 2008). The function of ELGN enzymes is inhibited in the absence of oxygen, and the consequent formation of the HIF complex promotes the transcription of several hundred hypoxia-response genes in order to reduce the effects of hypoxia. The expression of egln3, hbae, egln2, egln1_L, cox (subunit 4 isoform 1 and subunit 8) and ddit4 was consistently highest in salmon alevins raised under severe hypoxic conditions (30% O 2 ).
Similarly, the significant hypoxia-induced expression of three hemoglobin genes in European sea bass was associated with stimulation of egln3 (Cadiz et al. 2017). Work with alveolar type II epithelial cell cultures showed that the regulation of hemoglobin gene expression is likely to be indirectly regulated by upstream transcriptional activation by the HIF complex (Grek et al. 2011). Salmonids possess multiple hemoglobin isoforms with distinct functional properties and several isoforms are produced at specific developmental stages (Andersen 2020). Embryonic hemoglobins are involved in oxygen transport during early developmental stages (Maruyama et al. 2002), and salmon hbae was one of eleven hemoglobin subunit transcripts that were differentially expressed across the three O 2 treatments. Moreover, differential regulation of hemoglobin isoforms in response to hypoxia has been observed in various other fish species, such as red drum (Sciaenops ocellatus) and the Lake Victoria cichlid Haplochromis ishmaeli. This so-called 'hemoglobin isoform switching' is associated with increasing Hb-O 2 affinity and may safeguard oxygen transport in fish during hypoxic events (Rutjes et al. 2007;Pan et al. 2017).
Hypoxic regulation of the HIF axis has been shown to influence gene regulation and physiology of several other downstream biological pathways. Here we show that cytochrome c oxidase cox4-1 and cox8 together with ddit4 were significantly up-regulated in expression as oxygen treatment concentrations decreased (File S1). DDIT4 promotes survival under chronic stress conditions such as hypoxia, and ddit4 expression is induced by HIF-3 in human endothelial cells during prolonged severe hypoxia (Hay and Sonenberg 2004;Janaszak-Jasiecka et al. 2016). Episodic hypoxic treatment of European sea bass larvae was associated with altered metabolism in juvenile stage fish (Vanderplancke et al. 2015) and chronic exposure of zebrafish parents to hypoxia for several weeks led to significantly increased hypoxic resistance and body length in larval offspring (Ho & Burggren 2012). Hypoxic treatments of rainbow trout embryo also affect subsequent growth and glucose metabolism (Liu et al. 2017a, b). Exposure of Atlantic salmon embryos to hypoxia prior to hatching resulted in improved hypoxia tolerance of the hatchling, while hypoxia treatments during alevin development did not affect hypoxia tolerance of the alevin (Wood et al. 2019). Hence, the effects of hypoxia seem to depend on the frequency and duration and what developmental the treatment is applied and when the effects are measured. Further studies are needed to assess the hypoxic tolerance of treated fish at later life stages in response to sustained low oxygen treatment, and whether this might condition the fish for improved tolerance to low oxygen environments.

CONCLUSION
Our study has demonstrated that the use of hypoxia as a treatment during the early life stages of Atlantic salmon has a significant effect on DNA methylation and upregulation of the expression of hypoxiaresponse genes that play a key role in HIF-mediated oxygen homeostasis. These results, along with other salmonid studies reporting increased tolerance in hatched larvae after embryonic hypoxia treatment, demonstrates the potential for environmental stimuli to be utilized as treatments to shape gene regulation in farmed Atlantic salmon.
As hypoxia in the farming environment of Atlantic salmon is a key contributor to mortality (Randall et al. 1982;Iversen et al. 2005) continuing research to improve hypoxic tolerance is necessary. This and other studies have verified that low oxygen treatments can trigger an epigenetic response that upregulates genes in the HIF-1 hypoxia Figure 9 Egln3 expression compared to DNA methylation in the promoter. Relative egln3 RNA expression in salmon raised in 30% O 2 , 60% O 2 or 100% O 2 against DNA methylation (%) at the CpG sites 6 (A) and 8 (B) of the egln3 gene promoter in alevins at 700DD (n = 36). response pathway. Whether or not this has the potential to improve tolerance to hypoxia, and the effects of such treatments on the response of the fish to other environmental stressors (e.g., immune response to disease), general survival and other side-effects will need to be investigated more thoroughly.