Gene expression and DNA methylation changes in response to hypoxia in toxicant-adapted Atlantic killifish (Fundulus heteroclitus)

Adaptation to environmental stressors is critical for survival in rapidly changing ecosystems. Understanding the physiological and molecular responses that underlie adaptive mechanisms is essential for predicting organismal sensitivity. Fish populations, particularly those inhabiting coastal waters, often face multiple environmental challenges simultaneously, which can compound the stress on these organisms. One such common stressor is hypoxia – low oxygen levels in the environment – which is increasingly documented in coastal regions due to anthropogenic activities leading to excess nutrient loading, harmful algal blooms and climate change (Wallace and Gobler, 2021). The presence and accumulation of anthropogenic chemicals in coastal ecosystems poses an additional, co-occurring threat to the health of fish populations. The ability of fish to cope with hypoxia is mediated through a range of physiological, transcriptional, and epigenetic mechanisms. However, very little is known about how populations chronically exposed to toxicants respond to secondary stressors such as hypoxia.

The Atlantic killifish (Fundulus heteroclitus) is one of the most abundant estuarine fish distributed along the East coast of the USA, where they play an ecologically important role in the food webs (Able et al., 2007; Kimball and Able, 2007). Their ability to tolerate wide changes in environmental conditions, including temperature, salinity, oxygen and pH, have made them an ideal model species to investigate the biochemical, physiological and evolutionary basis of environmental adaptation (Borowiec et al., 2020, 2018; Borowiec and Scott, 2020). Some populations of killifish are also valuable models for understanding the mechanisms of evolved resistance to toxicants (Harishchandra et al., 2024; Nacci et al., 2010; Reid et al., 2016). Populations of killifish inhabiting contaminated coastal waters along the North Atlantic US coast have evolved resistance to some contaminants representing major categories of aryl hydrocarbon pollutants, such as polynuclear aromatic hydrocarbons (PAHs), and halogenated aromatic hydrocarbons such as polychlorinated biphenyls (PCBs), 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD, ‘dioxin’) and other dioxin-like compounds (DLCs) (Nacci et al., 2010). This evolved resistance involves alterations in signaling through the aryl hydrocarbon receptor (AHR), a ligand-activated transcription factor that forms a heterodimer with aryl hydrocarbon receptor translocator (ARNT), binding to dioxin-response elements to regulate the expression of target genes. While there is a great deal of understanding about the physiological and biochemical basis of adaptation to a variety of environmental conditions including toxicants in this species, very little is known about the impact of resistance to toxicants on their ability to respond to subsequent stressors such as hypoxia (Jayasundara et al., 2017; Jasperse et al., 2023; Lindberg et al., 2017).

The hypoxia-inducible factor (HIF) signaling pathway is a critical cellular response mechanism that enables organisms to respond to hypoxic conditions. Under normal oxygen levels (normoxia), HIF-α subunits (mainly HIF-1α and HIF-2α) are hydroxylated by prolyl hydroxylase enzymes, marking them for degradation via the von Hippel-Lindau (VHL) ubiquitin-proteasome pathway (Fong and Takeda, 2008; Kallio et al., 1999). This prevents the accumulation of HIF-α under normoxia. Under hypoxic conditions, the activity of prolyl hydroxylases is inhibited due to a lack of oxygen, leading to the stabilization of HIF-α (Ivan and Kaelin, 2017). Once stabilized, HIF-α translocates into the nucleus, where it dimerizes with ARNT, also known as HIF-1β. This HIF-α/ARNT complex binds to hypoxia response elements (HREs) and regulates the transcription of target genes (Kindrick and Mole, 2020; Li et al., 2023). Some of the target genes regulate processes such as angiogenesis (e.g. VEGF), erythropoiesis, glucose metabolism (e.g. GLUT1), and anaerobic metabolism (e.g. LDHA) (Luo et al., 2022). Similar responses were observed in killifish, suggesting conserved physiological and molecular mechanisms across species (Borowiec et al., 2020, 2018, 2015; Borowiec and Scott, 2020, 2021; Du et al., 2016; Lau et al., 2021; Ridgway and Scott, 2023).

The AHR and HIF pathways exhibit crosstalk primarily through their shared use of the heterodimerizing partner ARNT (Vorrink and Domann, 2014). Since ARNT is a limiting factor, competition between AHR and HIF for ARNT can influence the balance of responses to environmental toxins and hypoxic stress. This crosstalk may result in altered cellular outcomes, particularly in situations where both pathways are activated simultaneously, such as under environmental stress. One study tested this hypothesis in Atlantic killifish by exposing them to a dioxin-like PCB for 3 days, followed by a hypoxia challenge (Kraemer and Schulte, 2004). Prior PCB exposure disrupted the classical hypoxia response by increasing hepatic glycolytic enzyme activity, suggesting that dioxin-induced AHR activation could limit ARNT availability for the hypoxia response.

The objective of this study was to investigate the impact of evolved resistance to toxicants on response to acute hypoxia. We characterized the hepatic gene expression and DNA methylation patterns in response to two levels of hypoxia in two distinct populations of Atlantic killifish (Fundulus heteroclitus). One population is from Scorton Creek, Sandwich, MA, USA (SC), a relatively pristine environment, and is considered sensitive to environmental toxicants (Bello et al., 2001). The other population originates from New Bedford Harbor, MA, USA (NBH), a superfund site heavily contaminated with dioxin-like PCBs, where the killifish have evolved resistance to toxicants. The NBH population represents a unique case study in how toxicant-adapted organisms might exhibit trade-offs in their ability to cope with additional stressors, such as hypoxia. We hypothesized that fish from NBH would exhibit altered responses to acute hypoxia compared to the sensitive Scorton Creek (SC) population. Specifically, we predicted that the NBH fish would show differential hepatic gene expression and DNA methylation patterns in response to hypoxia, potentially compromising their ability to mount an optimal response to low oxygen conditions compared to the SC population.

To ensure that NBH and SC fish exhibit toxicant resistant and sensitive phenotypes, we exposed the embryos from both populations to PCB126 and quantified cyp1a gene expression using quantitative real-time PCR. SC embryos exposed to PCB126 induced cyp1a gene expression several hundred-fold compared to the DMSO group suggesting sensitivity to dioxin-like PCBs. In contrast, NBH fish showed a muted response to PCB126 exposure suggesting resistance to dioxin-like PCBs (Fig. 1B).

Neither mild nor severe hypoxia exposure had any effect on the loss of equilibrium during the 6-h exposure period in fish from both populations. Upon initial transfer into the hypoxia chamber, fish from both populations exhibited a rapid swimming response for the first 10-15 min. This was followed by a noticeable reduction in swimming activity accompanied by rapid ventilation (opercular movements). By the end of the exposure period, the fish were consistently found at the bottom of the container, exhibiting slow opercular movements but there was no loss of equilibrium as evidenced by their ability to maintain position and coordinated movement.

Strand-specific RNA sequencing of NBH and SC samples yielded an average of 17.2 million reads per sample. Of these, 83% of the reads were uniquely mapped to the genome. The summary of mapping statistics and read counts for annotated genes is provided in the supplementary information (RNAseq_supplementaryInformation.xlsx). Principal component coordinate analysis revealed one NBH hypoxia sample to be an outlier, which was omitted from statistical analysis (Fig. S2).

We observed a dose-dependent effect of hypoxia on differential gene expression in SC fish. Exposure to mild hypoxia revealed 2241 differentially expressed genes (DEGs), with 1170 upregulated and 1071 downregulated. In response to severe hypoxia, 4191 DEGs were observed, with 2221 upregulated and 1970 downregulated. A total of 1794 DEGs were shared between the two hypoxia groups, with 980 upregulated and 814 downregulated genes (Fig. 2A). Gene Ontology (GO) analysis of DEGs from mild and severe hypoxia treatment groups revealed overrepresentation of GO molecular function (MF) terms related to ATPase activity, RNA binding, proteosome and extracellular matrix functions. The list of top 10 overrepresented GO:MF terms among up- and downregulated DEG in SC mild and severe hypoxia groups are shown in Fig. 3.

In NBH fish, exposure to mild and severe hypoxia elicited differential expression of 3328 and 2423 genes, respectively. Among the 3328 genes differentially expressed in response to mild hypoxia, 1717 were upregulated and 1611 genes were downregulated. Whereas in response to severe hypoxia, 1272 of the 2423 DEGs were upregulated, and 1151 genes were downregulated. Comparison of the upregulated DEGs from the mild and severe hypoxia groups revealed 847 genes shared between the two hypoxia groups (Fig. 2B), while a similar comparison of downregulated DEGs revealed 796 shared DEGs.

GO analysis of mild and severe hypoxia upregulated DEGs showed overrepresentation of GO:MF terms related to mRNA splicing, translation, and proteasomal degradation. The terms enriched among downregulated DEGs include specific pathways such as cell adhesion and extracellular matrix functions. The top 10 overrepresented terms among up and downregulated DEGs are shown in Fig. 4.

Comparison of NBH and SC control groups revealed differential expression of 307 genes. Among them 159 and 148 are up- and downregulated in NBH, respectively, in comparison to SC. The GO terms enriched among these genes are shown in Fig. S3. Comparison of mean expression (log counts per million) of all up- and downregulated DEGs pooled from both populations and hypoxia treatment showed a dose-dependent response in SC fish, whereas NBH fish have a more muted gene expression response to severe hypoxia (Fig. 5).

We compared hypoxia treatment groups between populations and identified 249 DEGs under mild hypoxia in NBH compared to SC. Of these, 127 genes were upregulated, and 122 were downregulated. Under severe hypoxia, there were 831 DEGs in NBH relative to SC, with 441 genes upregulated and 390 downregulated. We pooled the DEGs of mild and severe hypoxia for GO analysis mainly due to the small number DEGs between populations exposed to mild hypoxia. DEGs that were upregulated in NBH fish were enriched MF terms such as oxidoreductase activity, iron and heme binding, and downregulated genes were enriched in MF terms that included protein tyrosine kinase and phosphatase activity, nucleic acid and chromatin binding, and electron transport chain pathways (Fig. 6). The full list of DEGs and GO analysis results are available in the supplementary information. This analysis largely supports the findings observed within the population-level hypoxia treatment effects (described above), showing that NBH fish, on average, exhibit a muted transcriptional response to hypoxia compared to SC fish.

Reduced representation bisulfite sequencing yielded 547.1 million total paired reads (14.4 to 33.3 million paired reads per sample). Of these reads, 544.5 million (99.5%) remained after quality-trimming. From the trimmed reads, 88.6% to 94.5% total reads were aligned to the F. heteroclitus genome. The detailed list of mapping statistics, raw and trimmed reads per sample is provided in the supplementary information (RRBS_Supplementary_information.xlsx). A total of 439,469 CpGs (5.4% of the 8,094,243 CpGs in the F. heteroclitus genome) had between 10× and 500× coverage in at least one sample after CpG filtering. Among them, 148,752 CpG sites (35.2% methylated and 64.8% unmethylated) were present in all samples and were used in downstream analysis.

There was no statistically significant difference in global methylation level between NBH (28.2%) and SC liver samples (28.8%). Fig. 7 shows the mean CpG methylation density plots and total number of CpG sites in all the treatment groups. Only one DMR was identified between NBH and SC fish exposed to control conditions. This DMR is in chromosome 10 in an intron of SHISA6, a gene that encodes AMPA glutamate receptor subunit (Klassen et al., 2016). This DMR also overlaps with an annotated CpG island (identified as CpG island 77 in the genome) and is hypomethylated in NBH in comparison to SC fish.

Hypoxia exposure did not substantially alter global DNA methylation levels in both populations. In SC fish, global CpG DNA methylation levels were 28.8% in the controls and mild hypoxia, and 29.2% in severe hypoxia group. No DMR were observed in response to mild or severe hypoxia in this population. In NBH fish, global DNA methylation level was 27.9% in the control group, 27.2% in mild hypoxia, and 29.7% in severe hypoxia. Comparison of control and mild hypoxia groups revealed 10 DMR, all of which were hypermethylated. Similar comparison between control and severe hypoxia revealed 59 DMRs. Among them five were hypomethylated and 54 were hypermethylated (Fig. 8). The majority of the NBH DMR were found in CpG islands annotated in the genome. Only two DMR were shared between the mild and severe hypoxia groups, and they were hypermethylated. The genomic coordinates of these DMR are provided in supplementary information. We did not observe any significant correlation between DNA methylation in DMRs and the expression level of the associated genes (Figs S4 and S5)

The findings from this study demonstrate a shift in physiological responses in fish adapted to toxicants and how this shift could compromise their capacity to respond to secondary stressors such as hypoxia. Our results show that both toxicant-sensitive (SC) and resistant killifish (NBH) respond to hypoxia but there are considerable differences in their responses as measured by gene expression and DNA methylation patterns. Atlantic killifish from SC show dose-dependent changes in gene expression patterns in response to hypoxia and no changes in DNA methylation patterns, whereas the fish from NBH showed a muted gene expression response to severe hypoxia suggesting compromised ability to mount a stress response. However, NBH fish showed a modest but significant DNA methylation changes in response to hypoxia. Together, these results reveal different response mechanisms and strategies to cope with hypoxia between toxicant-sensitive and toxicant-resistant fish populations.

Transcriptional responses to hypoxia are well-documented in fish species (Qin et al., 2023; Shang et al., 2022a,b; Zhao et al., 2022; Feng et al., 2022; Beck et al., 2016), including Atlantic killifish (Flight et al., 2011; Townley et al., 2017; Thomas and Kinsey, 2024). Our findings demonstrate that both SC and NBH killifish exhibit responses to hypoxia, though their reactions differ significantly. The hypoxia-inducible transcription factors (HIF-1α, HIF-2α, HIF-3α) play a major role in the transcriptional activation of the hypoxic response (Lisy and Peet, 2008). The expression of these genes under hypoxic conditions depends on the intensity and duration of hypoxia exposure. Interestingly, we did not observe any changes in the expression of HIF1α and HIF2α genes in response to mild or severe hypoxia in either of the populations. However, we observed an increased expression of HIF3α only in NBH fish (mild hypoxia −log FC 1.75, FDR 7.45E-04; severe hypoxia −log FC 1.54, FDR 4.07E-03). A previous study in a closely related species, Fundulus grandis, showed no significant differences in mRNA expression of any of these HIF genes in response to hypoxia (1 mg oxygen l⁻¹) after 6 and 24 h post-exposure (Murphy et al., 2023). However, studies in mammalian cell culture systems have shown differences in temporal profiles of HIF genes in response to hypoxia. In human endothelial cells, HIF-1α expression is maximal after 4 h of hypoxia and then is dramatically reduced by 8 h (Jaskiewicz et al., 2022a,b). In contrast, HIF-2α is maximal at 8 h and remains elevated up to 24 h (Jaskiewicz et al., 2022a,b; Moszynska et al., 2022). Unlike HIF-1α and -2α, moderate hypoxia exposure induced HIF3α mRNA expression within 2 h (Jaskiewicz et al., 2022a,b). However, very little is known about the functional role of HIF3α in hypoxia responses. The upregulation of HIF3α in response to hypoxia only in NBH fish suggests that adaptation to a contaminated environment may have altered its regulation, with potential implications for the hypoxia response. Further studies are needed to understand the roles of all three HIFs under different hypoxic conditions and at different life stages.

Another set of HIF pathway genes that were differentially expressed are the prolyl hydroxylase domain-containing proteins (PHD) (Fong and Takeda, 2008). There are three 2-oxoglutarate-dependent PHD proteins (PHD1, PHD2, and PHD3, encoded by EGLN2, EGLN1, and EGLN3, respectively) that are involved in HIFα ubiquitination and proteasomal degradation under normoxic conditions (Hoffman et al., 2001). Hypoxia has been shown to induce the expression of EGLN1 and -3 mRNAs, but not EGLN2, in several cell types (Ivan and Kaelin, 2017). In most cases, the induction of EGLN3 mRNA is much more prominent than that of EGLN1. We observed a several-fold increase in the egln3 expression (SC: log FC 8.77, FDR 1.82E-06; NBH: log FC 3.51, FDR 0.028) and modest but significant elgn1a upregulation in response to severe hypoxia in both populations (SC: log FC 1.33, FDR 9.59E-04; NBH: log FC 0.9, FDR 0.03). We also observed upregulation of egln2 (logFC 1.33, FDR 9.59E-04) in response to severe hypoxia in SC fish, while it was not significant in NBH fish. Although very little is known about the role of egln2 in oxygen sensing and hypoxia tolerance, the differences in expression between the two populations are intriguing. It has been suggested that increase in EGLN activity during hypoxia acts as a regulatory feedback loop for fast elimination of HIF after reoxygenation (del Peso et al., 2003).

In addition to the direct regulation of HIF transcription by hypoxia, multiple signaling pathways have been shown to play a role in the regulation of HIF gene expression in a variety of model systems. These include the PI3K-mTOR, interleukin-6 (IL-6), ERK, and MAPK signaling pathways (Luo et al., 2022). We observed differential expression of several genes associated with these signaling pathways only in NBH fish, suggesting a role in adaptation to toxicants. Notable among them include serine/threonine-protein kinase mTOR (mtor), phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit alpha (pik3r1) and RAC-beta serine/threonine-protein kinase (akt2). Despite the lack of an increase in HIF-1α expression, the differential expression of these signaling pathways suggests that NBH fish use crosstalk between multiple signaling pathways to cope with hypoxia in a way that appears to be distinct from SC fish. There is growing evidence that HIF-, PI3K-mTOR and ERK-MAPK signaling pathways act in an integrated way, influencing downstream pathways that affect cellular metabolism and survival during hypoxia (Luo et al., 2022). The differential expression of genes involved in the mTOR signaling pathway suggests disruptions in protein folding and secretion, triggering the unfolded protein response (UPR), which aids in restoring endoplasmic reticulum homeostasis (Wouters and Koritzinsky, 2008). The role of UPR in oxidative stress and hypoxia is well-established (Guan et al., 2024) and toxicant-resistant NBH fish would be an ideal model to investigate the effects of chronic exposure to toxicants on protein endoplasmic reticulum (ER) stress.

While responses to hypoxia are well documented across various organisms, encompassing metabolic changes, angiogenesis, and vascularization (Borowiec et al., 2018, 2015; Borowiec and Scott, 2020; Luo et al., 2022), our results, based on a 6-h hypoxia exposure, revealed gene expression changes that were more focused on transcription, translation, and cell cycle-related genes, rather than the classical hypoxia responses. This contrasts with studies in in vitro mammalian systems, where hypoxia target genes are differentially expressed in shorter timescales (Wu and Yotnda, 2011). The lack of similar responses in our study could be due to differences in the timescales of transcription and translation between fish maintained at 20°C and mammalian cells held at 37°C (Wu and Yotnda, 2011; Wenger et al., 2015).

We analyzed population differences in response to hypoxia using two approaches: first, by comparing hypoxia responses within each population, and second, by comparing hypoxia treatments between populations. Both approaches yielded similar conclusions. SC fish exhibited dose-dependent changes in gene expression, indicating a gradual and proportional response to varying levels of hypoxia. In contrast, NBH fish demonstrated drastic gene expression changes in response to mild hypoxia, but a muted response to severe hypoxia. This pattern suggests that NBH fish have a lower tolerance threshold for hypoxic conditions and are unable to induce a robust transcriptional response to severe hypoxia. This could be potentially due to higher energetic costs for survival in contaminated environments, leading to compromised ability to respond to additional stressors (Jayasundara et al., 2017; Jasperse et al., 2023).

In SC fish, mild hypoxia exposure caused an overrepresentation of genes associated with efflux transporters, transcription factor activity, mRNA and nucleic acid binding, and proteasomal functions. All these functions have been previously shown to be altered by hypoxia (Kallio et al., 1999; Vadlapatla et al., 2013; Thews et al., 2008; Yfantis et al., 2023; Chee et al., 2019). For instance, efflux pumps belonging to the ATP-binding cassette (ABC) superfamily of membrane transporters are known to play a significant role in cellular protection against oxidative stress (Patak et al., 2011). Several genes within the ABC family of efflux transporters are expressed in the liver and are involved in the transport of glutathione and glucuronide conjugates (Wlcek and Stieger, 2014; Kipp and Arias, 2002; Chan and Vandeberg, 2012). These findings suggest that mild hypoxia increases reactive oxygen species (ROS), with glutathione playing a crucial role in neutralizing free radicals and protecting the liver. The glutathione conjugates are subsequently eliminated from hepatocytes by the efflux pumps, contributing to cellular detoxification (Chan and Vandeberg, 2012). Additionally, mild hypoxia caused the upregulation of genes associated with RNA polymerase II (RNAPII) general transcription factor activity and the downregulation of cis-regulatory sequence-specific DNA binding RNAPII activity. This is not surprising, given the evidence that hypoxia affects the transcription of hundreds of genes (Kindrick and Mole, 2020). Under normoxic conditions, almost all HIF target genes display an open chromatin structure and harbor transcriptionally active but paused RNAPII (Soliman et al., 2024). Changes in the expression of RNAPII activity genes in response to hypoxia suggests the release of paused polymerase II into productive RNA synthesis by recruiting various coactivators, repressors, and chromatin remodelers, resulting in either the activation or inhibition of transcription of target genes (Soliman et al., 2024).

In response to severe hypoxia, SC fish showed an overrepresentation of genes associated with RNA helicases. They play an important role in cellular RNA metabolism, including transcription, pre-mRNA splicing, RNA export, storage, decay, and translation. Recently, RNA helicases have been shown to be involved in many biological processes, including DNA damage repair, cellular stress response, hypoxia and antiviral defense (Wang et al., 2021; Cai et al., 2017). Another major overrepresented group of genes in response to severe hypoxia in SC fish are the RNA binding proteins (RBPs) – important players in mRNA turnover (decay and stabilization) and translation (Mitchell and Tollervey, 2001). RNA helicases and RBPs are essential for the adaptive cellular response to hypoxia (Wang et al., 2021; Masuda et al., 2009). While there is very little understanding of the role of RBPs in environmental model species such as killifish, the upregulation of these genes in response to hypoxia suggests highly conserved cellular mechanisms.

Interestingly, the genes and pathways that were upregulated in SC fish under severe hypoxia were similarly upregulated in NBH fish exposed to mild hypoxia. This suggests that the toxicant-adapted NBH fish are more sensitive to mild hypoxia, undergoing more drastic transcriptional changes. The upregulation of RNA helicases and RBPs indicates post-transcriptional regulation of pre-existing RNAs (Masuda et al., 2009), suggesting that NBH fish may reduce transcription under stress. Indeed, this was observed under severe hypoxia in NBH fish, where the number of DEGs was lower than it was in SC fish, indicating a drastic reduction in transcription.

Hypoxia exposure downregulated genes related to calcium signaling, oxidoreductase activity, and extracellular matrix (ECM) modeling proteins. These pathways were altered in both hypoxia treatments and across both populations, suggesting they are vital for cellular response to hypoxia. It is well established that hypoxia impairs mitochondrial respiration and ATP synthesis, leading to an increased production of reactive oxygen species (ROS) and calcium release from the endoplasmic reticulum into the cytosol (Luo et al., 2022). The resulting elevated cytosolic calcium levels cause increased calcium uptake into mitochondria and mitochondrial calcium overload, which in turn leads to mitochondrial depolarization and the initiation of cell death (Seta et al., 2004). The decreased expression of genes associated with calcium signaling suggests an adaptation toward hypoxia tolerance. Similarly, hypoxia involves a switch from oxidative phosphorylation to glycolysis, resulting in increased production of NADH and an imbalance in the NAD+ and NADH ratio, causing altered redox potential (Li et al., 2023; Janczy-Cempa et al., 2022). These conditions favor the overexpression of many redox enzymes, such as cytochrome P450 reductase and nitroreductases. Even though these genes were not altered, genes that are dependent on NAD+ or NADH and play critical roles in metabolism were downregulated. They include agmo (alkylglycerol monooxygenase), sc5d (sterol-C5-desaturase), hsd11b2 (11-β-hydroxysteroid dehydrogenase 2), bdh1 (d-beta-hydroxybutyrate dehydrogenase), and foxred1 (FAD-dependent oxidoreductase domain-containing 1) (Li et al., 2023). It remains to be determined whether the downregulation of these genes under hypoxia is adaptive or maladaptive.

Another significant group of genes that are commonly downregulated in response to hypoxia are ECM components genes. The ECM is a complex network of proteins and other molecules that provide structural and biochemical support to surrounding cells (Naba et al., 2012). Its composition and function are crucial for tissue integrity, cell behavior, and overall organism health (Gilkes et al., 2014). Hypoxia has been shown to cause significant effects on the ECM in various aquatic species (Mu et al., 2020; Rojas et al., 2024; Shang et al., 2022a,b; Chen et al., 2021). In addition, there is growing evidence suggesting the role of hypoxia and HIFs in reprogramming cancer cells by regulating extracellular matrix (ECM) deposition, remodeling and degradation, thereby promoting cancer metastasis (Bertout et al., 2008). We observed downregulation genes associated with collagen synthesis, which maintain the structural integrity of tissues. This could lead to weakened tissue architecture or cause tissue remodeling such as changes in vascularization, which can be adaptive to survive under hypoxic conditions. The molecular mechanisms associated with these changes could be either by direct regulation of ECM pathway genes by HIF proteins or indirectly by the induction of oxidative stress (Gilkes et al., 2014; Vujić et al., 2022). Overall, the effects of hypoxia on the extracellular ECM in killifish highlight the importance of this adaptive function in coastal fish species, which frequently encounter periodic hypoxic conditions.

ARNT serves as a shared dimerization partner for both the AhR and HIF-1α signaling pathways, raising the possibility of competition between these pathways for ARNT. This hypothesis has been explored primarily using in vitro models, which demonstrate that ARNT can act as a limiting factor (Vorrink and Domann, 2014; Chan et al., 1999; Zhang et al., 2022). For instance, activation of the hypoxia pathway has been shown to inhibit AHR agonist-induced cyp1a upregulation, while activation of the AHR pathway has been shown to have either an additive effect (Chan et al., 1999) or no inhibitory effect on the expression of hypoxia responsive genes (Zhang et al., 2022). Similar patterns have been observed in other systems (Nie et al., 2001).

In this study, NBH fish exhibit constitutively high cyp1a expression compared to SC fish (log FC 1.306, FDR 0.009), suggesting that parental exposure to dioxin-like PCBs had a cross-generational effect on gene expression patterns in laboratory reared offspring. If ARNT is indeed limiting, NBH fish would be expected to show a reduced capacity for hypoxia responses. Consistent with this, we observed a muted response to severe hypoxia in NBH fish, suggesting that ARNT availability may be constrained. However, further studies are required to elucidate AHR:ARNT and HIF:ARNT protein interactions in these populations to confirm this hypothesis. Basal expression levels of ARNT genes (arnt and arnt2) do not differ significantly between the two populations, though arnt is differentially expressed in response to hypoxia in both. These findings suggest that assessing ARNT interactions at the protein level will provide stronger evidence for its role as a limiting factor and its influence on environmental adaptation in NBH fish.

Epigenetic effects, particularly DNA methylation changes in response to hypoxia, are well documented across a variety of fish species (Ragsdale et al., 2022; Thoral et al., 2022; Gu et al., 2024; Jones and Griffitt, 2022; Kelly et al., 2020). Hypoxia has been shown to cause both global and gene-specific alterations in DNA methylation patterns, which can influence gene expression, development, and stress response pathways. For instance, in species like zebrafish (Danio rerio) and Atlantic salmon (Salmo salar), hypoxic conditions have been associated with both hypermethylation and hypomethylation of key regulatory genes involved in metabolic adaptation and oxidative stress responses (Ragsdale et al., 2022; Kelly et al., 2020). These epigenetic changes are hypothesized to help fish cope with hypoxia by modulating critical pathways for survival, growth, and development.

To our knowledge, this is the first study where genome-wide DNA methylation patterns are profiled in Atlantic killifish. In this study, we hypothesized that resistant and sensitive populations of Atlantic killifish would exhibit distinct hepatic DNA methylation patterns, possibly reflecting their differential tolerance to dioxin-like PCBs. However, we did not observe distinct population differences in DNA methylation (NBH versus SC control groups), suggesting that parental exposure to dioxin-like PCBs did not have a detectable effect on DNA methylation in the offspring held under normoxic conditions.

Hypoxia exposure induced DNA methylation changes that were markedly different in the two populations. In the sensitive SC population, no differentially methylated regions (DMRs) were identified in response to either mild or severe hypoxia. In contrast, hypoxia resulted in significant changes in DNA methylation in the toxicant-resistant NBH population. Most of these DMRs are located in intronic regions and were enriched within CpG islands, which are known to play a critical role in regulating gene expression (Deaton and Bird, 2011). These findings suggest that the NBH population may have a more plastic epigenetic response to hypoxia compared to the SC population, which could reflect underlying differences in their adaptive capacities. While further research is needed to elucidate the functional consequences of these methylation changes, this study provides new insights into the epigenetic mechanisms by which fish populations respond to environmental stressors like hypoxia.

The majority of hypoxia induced DMRs are hypermethylated, supporting previous observations that hypoxia cause DNA hypermethylation (Thienpont et al., 2016; Rosa Ng and Jain, 2016; Robinson et al., 2012). Associating DMRs with genes revealed that some of the DMRs are related to genes involved in sialylation, vascularization and development. Four DMRs were associated with the gene St6galnac3, which is involved in the sialylation pathway. Sialylation refers to the addition of sialic acid units to oligosaccharides and glycoproteins (Li and Ding, 2019). Sialic acids moieties act as bridging molecules, facilitating communication between cells and the extracellular matrix. Additionally, one DMR was associated with beta-chimaerin, a protein linked to vascularization (Ben Dhaou et al., 2022). Hypoxia has been shown to influence both sialylation and vascularization processes, particularly in cancer models (Ben Dhaou et al., 2022; Jones et al., 2018), suggesting that similar mechanisms may be involved in the hypoxic responses observed in this study. However, we did not observe any significant correlation between these DNA methylation changes and expression of the associated genes. None of the 58 genes associated with DMRs in NBH fish were differentially expressed suggesting a temporal lag in DNA methylation changes and gene expression. In addition, several studies have shown that the majority of the DNA methylation changes are not correlated with gene expression (Aluru et al., 2018; Lindner et al., 2021; Bogan and Yi, 2024). This suggests that epigenetic regulation of gene expression is multilayered, with many levels of control, involving DNA methylation, histone modifications and chromatin organization.

While differential methylation did not correlate with gene expression changes in NBH, we observed differential expression of several chromatin modifier genes, particularly histone lysine demethylases (KDMs), in response to hypoxia in both fish populations. This is expected, as KDMs belong to the family of 2-oxoglutarate-dependent dioxygenases, which function as oxygen sensors (Kindrick and Mole, 2020). The differentially expressed KDM genes include kdm1aa, kdm2aa, kdm2ab, kdm2ba, kdm3b, kdm4b, kdm5ba, kdm5bb, kdm5c, kdm6a, kdm6ba, and kdm7aa. These genes have been shown to be directly regulated by hypoxia-inducible factors (HIFs), linking their expression to the cellular response to hypoxic stress (Salminen et al., 2016; Lee et al., 2014). Given the differential expression of histone lysine demethylases in response to hypoxia, it would be intriguing to investigate the role of chromatin modifiers in the adaptation to hypoxia in environmental species, as they may be key regulators of gene expression in low-oxygen environments.

This study highlights the complex and distinct physiological and epigenetic responses of Atlantic killifish populations adapted to toxicants when exposed to hypoxia. Our findings suggest that the capacity to respond to secondary stressors such as hypoxia may be altered in populations adapted to environmental contaminants. As expected, the toxicant-sensitive SC fish displayed a dose-dependent response to hypoxia exposure. However, the toxicant-resistant NBH fish, exhibited muted transcriptional responses but a more pronounced DNA methylation response to severe hypoxia, suggesting different molecular mechanisms in this population. Importantly, the differential DNA methylation patterns in response to hypoxia between the two populations indicate differences in epigenetic plasticity, which needs further investigation. A limitation of this study is the focus on liver tissue, which may not fully represent systemic responses to hypoxia. Additionally, physiological differences between the NBH and SC populations, shaped by their distinct habitats, may confound the interpretation of molecular responses. Overall, this research provides valuable insights into the diverse molecular mechanisms by which fish populations, with different environmental histories, respond to hypoxia and highlights the need for further exploration of epigenetic and chromatin-level responses in the context of environmental adaptation.

The animal husbandry and experimental procedures used in this study were approved by the Animal Care and Use Committee of the Woods Hole Oceanographic Institution. Mature adult male and female killifish from Scorton Creek (SC; Sandwich, MA, USA) and New Bedford Harbor (NBH; New Bedford, MA, USA) were collected during low tide in the month of June 2018 using minnow traps, as described previously (Karchner et al., 1999). The water temperature and salinity at the time of collection in NBH was 27°C and 34 parts per thousand (ppt) and in SC was 25°C and 33 ppt respectively. In NBH, fish were collected from the banks of the Acushnet River estuary (latitude 41.670, longitude −70.915) during receding tide. Fish from SC were collected from tidal pools (latitude 41.746, longitude −70.428). Fish were brought to the lab and maintained in the Redfield Laboratory (WHOI) with continuous flow-through seawater (SW) at 18–20°C, saturated dissolved oxygen (21% oxygen saturation or 7.21 mg oxygen l⁻¹), and 14 h:10 h light:dark photoperiod conditions.

F1 generation of embryos from SC and NBH were obtained by in vitro fertilization following established protocols (Karchner et al., 1999). Briefly, four or five female fish from SC or NBH (15.5±1.2 g mean wet mass) were lightly anesthetized with Tricaine (MS222; buffered with sodium bicarbonate, Sigma-Aldrich, St. Louis, MO, USA) and oocytes were obtained for in vitro fertilization by gently squeezing the abdomen. Oocytes were collected in glass Petri dishes with filtered SW (30 parts per thousand; ppt). Milt was obtained by euthanizing 2-3 mature males (13.5±1.1 g mean wet mass) from the same population in MS222, dissecting out the gonads, and chopping them with a scalpel blade in seawater. A few drops of milt were added to the oocytes for fertilization. Approximately 20 min after the addition of milt, embryos were rinsed with filtered SW to remove any excess sperm. Fertilized embryos were reared at 23°C under 14 h:10 h light:dark photoperiod conditions until hatching. Larvae were raised in 2-gallon aquarium tanks in aerated seawater for 6 months. During larval rearing, fish were fed brine shrimp daily and water was exchanged every 2-3 weeks. Oxygen concentration was measured in the tanks once every 2-3 days and oxygen saturation was above 20% throughout the rearing period.

A subset of F1 generation embryos from NBH and SC were exposed to either DMSO (0.01%; vehicle control) or 50 nM PCB126 from 4 h post-fertilization to 5 days post-fertilization through water borne exposure. Exposures were conducted in glass Petri dishes in 50 ml of seawater at 23°C under 14:10 h light:dark photoperiod regimen. Each treatment condition consisted of 30 embryos. At the end of the exposure, pools of six embryos were flash frozen in liquid nitrogen for RNA isolation (five replicates per treatment). Total RNA was isolated using the Aurum total RNA fatty and fibrous tissue kit (BioRad, Herculus, CA, USA) with in column DNase treatment. cDNA was prepared with using iScript cDNA synthesis kit (BioRad, CA, USA) and used in qRT-PCR to quantify cyp1a expression. β-actin was used as a housekeeping gene. The qRT-PCR primers used to amplify cyp1a were 5′-CTTTCACAATCCCACACTGCTC-3′ (forward primer) and 5′-GGTCTTTCCAGAGCTCTGGG -3′ (reverse primer). β-actin primers used were 5′-TGGAGAAGAGCTACGAGCTCC-3′ (forward primer) and 5′-CCGCAGGACTCCATTCCGAG-3′ (reverse primer). The PCR conditions used were 95°C for 3 min and 95°C for 15 s/66°C (cyp1a and β-actin) for 1 min (40 cycles). At the end of each PCR run, a melt curve analysis was performed to ensure that only a single product was amplified. Three technical replicates were used for each sample. Relative expression was normalized to that of β-actin {2^−ΔCt; where ΔCt=[Ct(cyp1a)−Ct(β-actin)]}. Cyp1a mRNA expression levels in response to PCB126 exposure were compared to respective population controls using one-way ANOVA (GraphPad Prism version 5.3). A probability level of P<0.01 was considered statistically significant. The embryos from this clutch were reared for 6 months and used in hypoxia exposure experiment described below.

Six-month-old killifish juveniles from SC (149±43 mg mean wet mass) and NBH (145±39 mg mean wet mass) were exposed to either mild (10% oxygen saturation, 3.46 mg O₂ l⁻¹; n=5 per population) or severe hypoxia (5% oxygen saturation, 1.72 mg O₂ l⁻¹ ; n=5 per population) for 6 h. These two hypoxia levels were chosen based on preliminary experiments with the same cohort of fish, where the loss of equilibrium (LOE) was assessed under 1% and 5% oxygen saturation in both NBH and SC fish. Six hours of exposure to 5% oxygen saturation did not cause LOE in fish from either population (n=5 individual fish per population), whereas 1% oxygen saturation caused LOE within 6 h in fish from both populations. The details of the results from the preliminary experiments are provided in Fig. S1.

Hypoxia exposure set up included Pyrex glass dishes (270 ml volume) equipped with oxygen sensor spots (PreSens Precision Sensing GmbH, Germany) placed inside hypoxia chambers (Stemcell Technologies) with pre-mixed air set to 5 or 10% oxygen pumped into the chambers continuously. A control group (normoxia, 20.9% oxygen saturation; n=5 per population) was maintained on the benchtop (Fig. 1). Prior to introducing the fish to hypoxia, 250 ml filtered seawater was added to the Pyrex glass dishes and was allowed to equilibrate overnight to ensure that the water had reached the respective treatment conditions. Oxygen levels in individual beakers were checked prior to introducing the fish using a FireString oxygen sensor (PyroScience, Germany) and were found to be at the treatment conditions in each of the individual dishes. At the start of the experiment, individual fish were quickly introduced into the Pyrex dishes and chambers closed immediately. Fish were maintained at treatment conditions for 6 h. At the end of the exposure period, oxygen levels were measured. Fish were rapidly anesthetized using a lethal dose of MS222 buffered with sodium bicarbonate. Livers were dissected and stored at −80°C until further analysis. Liver tissue was selected to quantify hypoxia-induced transcriptional and DNA methylation responses due to its central role in metabolism. Hypoxia has been widely shown to significantly impact metabolic pathways, and in toxicant-resistant killifish, studying the response to hypoxia offers valuable insights into how adaptation to toxicants influences metabolic homeostasis under low-oxygen conditions.

Simultaneous isolation of genomic DNA and total RNA from liver tissues was performed using the ZR-Duet DNA/RNA Mini Prep kit (Zymo Research, CA, USA). RNA was treated with DNase during the isolation process. DNA and RNA were quantified using the Nanodrop Spectrophotometer. The quality of DNA and RNA was checked using the Agilent 4200 and 2200 Tape Station systems, respectively. The DNA and RNA integrity numbers of all samples were between 9 and 10.

Libraries were constructed using Illumina stranded library preparation kit following manufacturer's protocol. Single end 50 bp reads were sequenced using Illumina HiSeq2500 platform. RNA sequencing library construction and sequencing were done at the Tufts University core facility.

Library preparation was performed using the Premium RRBS kit (Diagenode). In brief, 100 ng DNA from each sample were enzymatically digested by the restriction enzyme MspI at 37°C for 12 h. Following ends preparation, a different set of adaptors was added to each sample and adaptor ligation was performed by the addition of ligase. Size selection of adaptor-ligated DNA fragments was performed by Agencourt AMPure XP beads (Beckman Coulter) and the DNA was eluted in Resuspension buffer (Diagenode, catalog number C02030033). Part of the eluted sample was subjected to qPCR using 2× KAPA HiFi HotStart ReadyMix (Kapa Biosystems) for quantification and subsequent pooling per nine samples. The pooling was performed according to two parameters: the Ct value and the adaptor ID of each sample. The pooling was followed by a cleanup with AMPure XP beads to reduce the volumes. Bisulfite treatment was performed, and bisulfite-converted DNA was eluted twice in BS Elution buffer (Diagenode, catalog number C02030033). Part of the bisulfite converted library was used in qPCR for the determination of the optimal cycle number for the enrichment PCR. 2× MethylTaq Plus Master Mix was used for the amplification PCR and a last cleanup with AMPure XP beads followed. PCR product was run on an 2% agarose gel to remove adaptor dimers. The quality of the final libraries was checked on an Agilent 2100 High Sensitivity DNA chip. The concentration was determined by performing qPCR on the samples using a dilution of PhiX index3 as standard. Paired end 50 bp reads were sequenced on an Illumina HiSeq4000 platform by a commercial facility (NXT-Dx, Ghent, Belgium).

The F. heteroclitus genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_011125445.2/) was used for all analyses. Genome feature information was pulled directly from the genome to generate gene, coding sequence (CDS), exon, and lncRNA genome feature tracks using Gnomom, RefSeq, cmsearch, and tRNAscan-SE annotations. Chromosome name and length information was also extracted from the genome to generate additional feature tracks using bedtools version 2.31.1 (Quinlan and Hall, 2010). A non-coding sequence track was created by using the complement of the CDS track (complementBed). Similarly, the intergenic genome feature track was created using the complement of the gene track. The intersection (intersectBed) between the non-coding sequence and gene tracks were used to create an intron track. All genome feature tracks are available in an Open Science Framework repository (doi.org/10.17605/OSF.IO/NZRA8).

Raw data files were assessed for quality using FastQC version 0.11.9 (Andrews, 2010) prior to preprocessing. Preprocessing was done by trimming the adaptor sequences using Trimmomatic (version 0.25) and removing any reads with low sequence quality (Phred score <20) (Bolger et al., 2014). Trimmed sequence reads were mapped to the F. heteroclitus genome using the STAR aligner version 2.6.1d (Dobin and Gingeras, 2015). The number of reads mapped to annotated regions of the genome was obtained using HTSeq-count version 0.11.1 (Anders et al., 2015). Statistical analysis was conducted using edgeR version 3.40.2, a Bioconductor package (Robinson et al., 2010). Transcripts from all samples were compiled into a DGEList, and lowly expressed transcripts were filtered out using the filterByExpr function. Sample ordination was visualized using multidimensional scaling analysis with the ape version 5.8 package (Oksanen et al., 2022), revealing an outlying sample that was removed from subsequent analysis (Fig. S2). We used the quasi-likelihood model in edgeR (glmQLFTest) to perform differential gene expression analysis. Only genes with false discovery rate (FDR) of <5% were considered to be differentially expressed. Raw data has been deposited in gene expression omnibus (accession number GSE278569).

Functional annotation of DEGs was done using gene ontology (GO) biological process (GO:BP) and molecular function (GO:MF) terms. Identification of overrepresented GO terms (P value <0.05) among sets of DEGs was done using the enricher function in ClusterProfiler version 4.6.2 (Wu et al., 2021). The background gene list included all expressed transcripts in the filtered DGEList. Similar GO:BP or GO:MF terms were clustered based on the frequency of shared genes using the R package rrvgo 1.16.0 (Sayols, 2023). Representative parent terms from each cluster were chosen based on the lowest P value.

The Bisulfite Analysis Toolkit (BAT) (Kretzmer et al., 2017) was used for RRBS analysis. Prior to analysis, raw data was quality trimmed with TrimGalore! version 0.6.6 (Martin, 2011). Trimming was performed on non-directional (--non_directional) paired-end reads (--paired). An additional 2 bp were trimmed from the 3′ end of the first read and 5′ end of the second read (--rrbs). Sequence quality was assessed with FastQC version 0.11.9 (Andrews, 2010) and MultiQC version 1.11 (Ewels et al., 2016) after trimming.

Trimmed paired reads were aligned to the genome using BAT_mapping module specifying non-directional input (-F 2). Mapping statistics and methylation calling was done using BAT_mapping_stat and BAT_calling modules, respectively. CpG methylation data was filtered to retain only a minimum 10 and maximum 500 reads per sample (--MD_min 10, –MD_max 500 --CG). The data were sorted using bedGraphs (sortBed version 2.29.1; Quinlan and Hall, 2010) and merged into treatment-specific groups (BAT_summarize). Within each group, one sample was allowed to have missing data for a CpG locus (--mis1 1, --mis2 1). If data were missing for more than one sample at a particular CpG, it was not included in the downstream analysis. Chromosome lengths were specified (--cs) for merging methylation information accurately. BAT_overview was used to obtain average methylation rate per sample in each group, hierarchical clustering of sample methylation rates, distribution of CpG methylation, comparison of methylation rate between groups for common loci, and differences in mean methylation rate between groups. Differentially methylated regions (DMR) — defined as regions with at least 10 CpGs, a minimum methylation rate difference of 0.1, and q-value <0.05 were identified for each comparison using BAT_DMRcalling module. The closest genome feature to each DMR was characterized using closestBed.