Transcriptomic and Methylomic Analyses Show Significant Shifts in Biosynthetic Processes and Reduced Intrapopulation Gene Expression Variance in PAH-Adapted Atlantic Killifish

Abstract

Environmental contaminants pose a significant selection pressure across taxa, potentiating evolved resistance to chemicals. However, rapid evolution may alter molecular and physiological homeostasis leading to trade-offs. To elucidate molecular underpinnings of evolved chemical resistance, we compared liver gene expression and methylation profiles in polycyclic aromatic hydrocarbon (PAH)-adapted Atlantic killifish (Fundulus heteroclitus) in the Republic site (RP), Elizabeth River, Virginia with PAH-sensitive Kings Creek (KC) fish. We found 1607 differentially expressed and 2252 alternatively spliced genes between RP and KC, with highly enriched genes involving lipid and amino acid metabolism, respectively. While 308 genes had differentially methylated regions, only 13 of these genes were differentially expressed. The aryl hydrocarbon receptor 2b gene (ahr2b) was differentially methylated and expressed, as well as alternatively spliced signifying its critical role in mediating PAH tolerance. Notably, the intrapopulation coefficient of variation (CoV) was lower in 82% of 17,566 expressed genes in RP fish compared to KC fish. Among other pathways, these genes with low CoV were highly enriched in bioenergetic processes inferring reduced metabolic physiological variation as a population in RP fish. Altered metabolic gene expression and overall reduced gene expression variance in RP fish warrant further studies on fitness trade-offs including altered susceptibility to other stressors associated with rapid adaptation to anthropogenic pressures.

Keywords: AHR pathway; Atlantic killifish; Elizabeth River; epigenetics; gene expression; pollution; polycyclic aromatic hydrocarbons; transcriptomics.

Figure 1.

Atlantic killifish collection sites for the study. The blue dot in the upper figure indicates the Kings Creek (KC, pollution-sensitive fish) site, and the red dot represents the Republic (RP, pollution-tolerant fish) site in the Elizabeth River, Virginia. Data are adapted from Clark et al. (2013).

Figure 2.

Genes associated with the reduced aerobic metabolic processes in PAH-adapted Republic fish (RP), compared to Kings Creek (KC) killifish. (a) Ratio of the coefficient of variation (CoV) for each individual gene expressed in killifish from KC and RP sites. Each dot represents a unique gene. Genes with a KC-CoV:RP-CoV ratio higher than one (on and above the solid dark line) are indicated in red dots. These are genes with high intrapopulation variance in KC. Genes with a KC-CoV:RP-CoV ratio lower than one are shown in green color. These are genes with high intrapopulation variance in RP. (b) Biological processes and cellular components of GO enrichment terms of genes with a CoV ratio higher than 2 (genes shown in red dots, above the dashed dark line) in the KC fish population compared to RP fish.

Figure 3.

Summary of the DEGs (717 downregulated and 514 upregulated genes in RP fish compared to KC fish) and alternatively spliced genes in Atlantic killifish from the RP site compared to fish from the KC site. (a) Among the listed AHR paralogues, only ahr2b showed higher expression and downregulation in the RP population. (b) Biological processes and molecular functions of enriched GO terms of DEGs (FDR < 0.05, n = 1607) in RP fish compared to KC fish. (c) Network map of the functionally connected (protein−protein interaction) genes among the DEGs (msmo1: methylsterol monooxygenase 1, mvd: mevalonate diphosphate decarboxylase a, nsdhl: NAD(P)-dependent steroid dehydrogenase-like, tm7sf 2: transmembrane 7 superfamily member 2, dhcr24: 24-dehydrocholesterol reductase, and cyp51: cytochrome P450, family 51). (d) Cellular components, biological processes, and molecular functions of GO enrichment terms of shared differentially expressed and alternatively spliced genes. RP, Republic site; KC, Kings Creek site.

Figure 3.

Summary of the DEGs (717 downregulated and 514 upregulated genes in RP fish compared to KC fish) and alternatively spliced genes in Atlantic killifish from the RP site compared to fish from the KC site. (a) Among the listed AHR paralogues, only ahr2b showed higher expression and downregulation in the RP population. (b) Biological processes and molecular functions of enriched GO terms of DEGs (FDR < 0.05, n = 1607) in RP fish compared to KC fish. (c) Network map of the functionally connected (protein−protein interaction) genes among the DEGs (msmo1: methylsterol monooxygenase 1, mvd: mevalonate diphosphate decarboxylase a, nsdhl: NAD(P)-dependent steroid dehydrogenase-like, tm7sf 2: transmembrane 7 superfamily member 2, dhcr24: 24-dehydrocholesterol reductase, and cyp51: cytochrome P450, family 51). (d) Cellular components, biological processes, and molecular functions of GO enrichment terms of shared differentially expressed and alternatively spliced genes. RP, Republic site; KC, Kings Creek site.

Figure 4.

Summary of the weighted gene coexpression analysis (WGCNA) conducted for Atlantic killifish genes expressed in both RP and KC populations. (a) Distribution of killifish genes among each highly correlated gene module (module eigengenes, MEs). (b) Functional level (protein−protein interaction) map of the most highly correlated (module members (MM), genes in a module eigengene with absolute correlation >0.8) genes in the ninth module eigengene (msmo1: methylsterol monooxygenase 1, dhcr24: 24-dehydrocholesterol reductase, idi1: isopentenyl-diphosphate delta isomerase 1, sc5d: sterol-C5-desaturase, srebf 2: sterol regulatory element binding transcription factor 2, dhcr7: 7-dehydrocholesterol reductase, and cyp51: cytochrome P450, family 51) (note that all these genes were differentially downregulated in the RP population). (c) GO enrichment terms of module members of ninth module eigengene.

Figure 4.

Summary of the weighted gene coexpression analysis (WGCNA) conducted for Atlantic killifish genes expressed in both RP and KC populations. (a) Distribution of killifish genes among each highly correlated gene module (module eigengenes, MEs). (b) Functional level (protein−protein interaction) map of the most highly correlated (module members (MM), genes in a module eigengene with absolute correlation >0.8) genes in the ninth module eigengene (msmo1: methylsterol monooxygenase 1, dhcr24: 24-dehydrocholesterol reductase, idi1: isopentenyl-diphosphate delta isomerase 1, sc5d: sterol-C5-desaturase, srebf 2: sterol regulatory element binding transcription factor 2, dhcr7: 7-dehydrocholesterol reductase, and cyp51: cytochrome P450, family 51) (note that all these genes were differentially downregulated in the RP population). (c) GO enrichment terms of module members of ninth module eigengene.

Figure 5.

DEGs-DMGs integrated analysis. (a) Log fold change of DEGs (blue stems) and DMGs (brown and green bars). Fourteen genes were both DEGs (FDR < 0.05 and absolute log fold change >1) and DMGs (FDR < 0.05 and absolute methylation difference >25%). The directionality of DMGs is indicated in the thick bars, where hypermethylation is indicated in brown color bars and hypomethylation is indicated in green bars. Blue stems to the right of the dark vertical line indicate upregulation of a gene and to the left indicate downregulation. Note that the log fold change values were multiplied by 10 for easy visualization. (b) Intron hypermethylation in two key AHR paralogues, ahr1b and ahr2b. This is a graphical illustration based on actual gene region data extracted from the Atlantic killifish genome version 3.0.2 available at the Ensembl genome browser. (wdr47—WD repeat domain 47; vwc2—von Willebrand factor C domain containing 2, ugdh—UDP-glucose 6-dehydrogenase; phactr—phactr1—phosphatase and actin regulator 1; pax7b—paired box protein 7b; limd1—LIM domain containing 1; kdm6bb—lysine demethylase 6b; eml1—echinoderm microtubule associated protein like 1; dclk1a—doublecortin like kinase 1; ahr2b—aryl hydrocarbon receptor 2b)