Epigenetic inheritance of DNA methylation changes in fish living in hydrogen sulfide-rich springs.

Environmental factors can promote phenotypic variation through alterations in the epigenome and facilitate adaptation of an organism to the environment. Although hydrogen sulfide is toxic to most organisms, the fish Poecilia mexicana has adapted to survive in environments with high levels that exceed toxicity thresholds by orders of magnitude. Epigenetic changes in response to this environmental stressor were examined by assessing DNA methylation alterations in red blood cells, which are nucleated in fish. Males and females were sampled from sulfidic and nonsulfidic natural environments; individuals were also propagated for two generations in a nonsulfidic laboratory environment. We compared epimutations between the sexes as well as field and laboratory populations. For both the wild-caught (F0) and the laboratory-reared (F2) fish, comparing the sulfidic and nonsulfidic populations revealed evidence for significant differential DNA methylation regions (DMRs). More importantly, there was over 80% overlap in DMRs across generations, suggesting that the DMRs have stable generational inheritance in the absence of the sulfidic environment. This is an example of epigenetic generational stability after the removal of an environmental stressor. The DMR-associated genes were related to sulfur toxicity and metabolic processes. These findings suggest that adaptation of P. mexicana to sulfidic environments in southern Mexico may, in part, be promoted through epigenetic DNA methylation alterations that become stable and are inherited by subsequent generations independent of the environment.

Toxicants are present in a wide variety of ecological contexts including natural and man-made, and environmental toxicants are becoming increasingly abundant due to anthropogenic activities. While many chemicals are categorized by their mutagenic effects, it has become increasingly clear that indirect modes of action are as important for short-term (e.g., phenotypic plasticity) and long-term (e.g., adaptive) effects (1–4). Responses to the toxicity of environmental chemicals have been shown, in part, to be mediated by epigenetics (4–6). For example, exposure to heavy metals, including cadmium and arsenic, has been associated with changes in DNA methylation in mammals (7–10). However, the findings are often complex and contradictory. In vitro cadmium exposure experiments in rat liver cells suggested initial inductions of DNA hypomethylation occurred, while prolonged exposure led to DNA hypermethylation (7). A wide variety of environmental compounds, such as the agricultural fungicide vinclozolin and the herbicide glyphosate, have been shown to modify the location and abundance of DNA methylation in rats (11–13). Moreover, vinclozolin and glyphosate promote the epigenetic transgenerational inheritance of DNA methylation states and disease susceptibility in the third generation following transient exposures (6, 11–13). The ability of these environmental toxicants to shape phenotypic variation through epigenetic changes that are stable across generations has potential evolutionary consequences that remain largely unexplored.

Natural systems with environmentally derived toxicants provide a framework in which it is possible to study the evolutionary effects of toxicants. By comparing populations adapted to a toxicant to nonadapted ancestral populations, we can determine the effect of the toxicant on the evolution of the populations. To make predictions about the outcomes from environmental exposure, it is necessary to have a model stressor with clearly defined and predictable effects. Hydrogen sulfide (H2S) is an ideal toxicant to study the role of epigenetics in responding to environmental toxicants, because H2S exposure has clear effects on sulfide processing and energy metabolism (14).

H2S is one of the most toxic inorganic gases for metazoan organisms. It occurs both naturally (e.g., in deep-sea hydrothermal vents, marine sediments, and sulfide springs) and as a by-product of pollution and industrial processes [e.g., in habitats impacted by farming, tanning, paper manufacturing, sewage treatment, oil refining, and gas exploration and refining (15, 16)]. For humans, with a recommended industrial daily exposure of 10 ppm (293 µM) (National Institute for Occupational Safety and Health, 2020), H2S can have detrimental health effects, including headache and eye irritation, even at low concentrations of 2.5 to 5 ppm (73 to 146 µM) (16), and exposure to high concentrations (>1,000 ppm/29,343 µM) can lead to instantaneous death (17). The primary reason for H2S’s high toxicity is that it directly inhibits cytochrome c oxidase, which is Complex IV of the mitochondrial electron transport chain (16, 18). The inhibition of cytochrome c oxidase by H2S results in a shutdown of the electron transport chain and cellular ATP generation (16, 18).

While H2S is toxic at high concentrations, H2S is produced at low concentrations endogenously as a product of cysteine catabolism and by intestinal bacteria (19, 20), and physiologically relevant concentrations are likely in the nanomolar range (21, 22). Due to the endogenous production of H2S, most organisms are able to detoxify low concentrations of H2S, and there are known sulfide detoxification enzymes present in most metazoans (23). Sulfide oxidation to thiosulfate is mediated by sulfide quinone oxidoreductase, a dioxygenase, and sulfur transferase (24). H2S is highly membrane permeable (25) and plays an important role as a cell-signaling molecule in the cardiovascular and nervous systems (26). For example, H2S is involved in the regulation of vasodilation and inflammation (27). The effect of H2S on epigenetic changes has not been investigated.

Very few organisms are able to tolerate exposure to high H2S concentrations (16, 28). However, some fish in the livebearing family Poeciliidae live in habitats with naturally occurring sustained and high concentrations of H2S (29–34). The current study focuses on the species Poecilia mexicana, which inhabits H2S-rich springs in southern Mexico (Fig. 1) (35). H2S in these springs ranges from 30 µM to over 1,000 µM, depending on the site (35). The concentrations fluctuate little across seasons and years (30, 36). There are also closely related populations of P. mexicana residing in nonsulfidic habitats adjacent to sulfidic springs. The populations are genetically differentiated but considered the same species (37–39). Wild-caught fish from sulfide spring populations exhibit significantly higher sulfide tolerance than fish from adjacent nonsulfidic populations (35). We have previously shown that P. mexicana from sulfidic populations can survive in high levels of H2S by constitutively expressing high levels of important H2S detoxification genes (40–42), and some sulfidic populations have evolved a resistant cytochrome c oxidase (43). Moreover, there is heritable variation in gene expression, especially in key genes related to H2S toxicity and detoxification (44), some of which may be driven by expression differences of relevant microRNAs (45). Despite these observations, a key gap in our knowledge is the effect of H2S exposure on epigenetics and whether expression changes due to H2S exposure and local adaptation are potentially mediated by DNA methylation changes. Fish are known to have an increased number of DNA methylation enzymes, DNA methyltransferases (DNMTs), and respond to environmental insults through DNA methylation (46).

Fig. 1.

Collection site, P. mexicana fish morphology, and experimental design. (A) Collection sites in southern Mexico. Map from ref. 79, which is licensed under CC BY 3.0. (B) Images of P. mexicana adult male and female fish from sulfidic and nonsulfidic environments. (C) Experimental design.

The current study was designed to understand the long-term epigenetic changes that can occur in response to H2S adaptation. This study also sets the stage for future studies examining how short- and long-term environmental exposure may shape the epigenome and the role of epigenetic transgenerational inheritance in adaptation to H2S. The primary research objectives were to test whether 1) the abundance and distribution of differential DNA methylation regions (DMRs) differ between sulfidic and nonsulfidic populations for both sexes, when samples were derived from wild versus laboratory-reared individuals, and 2) the methylation differences between populations are consistent and related to living in H2S. We hypothesized that long-term H2S exposure in P. mexicana involved alterations in epigenetics (DNA methylation) and that those changes are stable over generations.

Results

Red blood cells (RBC) from male and female P. mexicana from H2S-rich habitats and nonsulfidic habitats in southern Mexico were collected to assess the environmental impacts of the sulfidic habitats on the fish epigenome (Fig. 1). The RBC in fish species contain nuclei, which allows DNA extraction from an easily purified single-cell type. Another set of P. mexicana from sulfidic and nonsulfidic habitats were collected as adults (pregnant females) in May 2013 and transported to a nonsulfidic laboratory environment and propagated for two generations in the nonsulfidic environment (Fig. 1). Adult F2 generation laboratory fish harvested in August 2015 were used for analysis and comparisons (Fig. 1). Genomic DNA from the nucleated RBC was isolated and used to identify differential DMRs between the sulfidic and nonsulfidic populations for males and females, separately. This experimental design allowed us to determine population differences in DMRs of wild fish populations and the stability of the epigenetic alterations in the laboratory propagated fish.

The RBC DNA was fragmented and the methylated DNA immunoprecipitated (MeDIP) with a methyl cytosine antibody to identify DMRs between the sulfidic and nonsulfidic populations. On average, similar numbers of DMRs were observed for the various comparisons between sulfidic and nonsulfidic populations (Table 1). Two-thirds of DMRs had one significant 100-base pair (bp) window, and the other third of DMRs exhibited multiple, significant 100-bp regions within the DMR (Table 2). The DMRs for each comparison are presented in . For both wild male and female comparisons between sulfidic and nonsulfidic populations, ∼50% of the DMRs had an increase in DNA methylation in the sulfidic population, and the remainder of the DMRs had a decrease in methylation (). Although a limitation of the study is the small number of fish analyzed, the robust and significant epigenetic data obtained support the observations presented.

Table 1.

DMR identification

Comparison	All window	Multiple window
Sulfidic Male Wild (SMW) vs. Nonsulfidic Male Wild (NMW)	1,049	440
Sulfidic Female Wild (SFW) vs. Nonsulfidic Female Wild (NFW)	1,461	606
Sulfidic Male Laboratory (SML) vs. Nonsulfidic Male Laboratory (NML)	1,619	641
Sulfidic Female Laboratory (SFL) vs. Nonsulfidic Female Laboratory (NFL)	1,451	544

The number of DMRs found in each comparison using P value threshold of P < 10−7. The All window (100 bp) column shows all DMRs. The Multiple window column shows the number of DMRs with multiple significant 100-bp windows within the DMR.

Table 2.

Number of DMRs with a given number of significant windows at p-value < 10−7

Comparison	Number of significant windows
	1	2	3	4	5	6	7	8	9	10	11	≥12
Sulfidic male wild (SMW) vs. nonsulfidic male wild (NMW)	609	240	76	45	23	17	10	4	8	5	3	9
Sulfidic female wild (SFW) vs. nonsulfidic female wild (NFW)	855	331	119	61	37	17	12	3	9	2	4	10
Sulfidic male laboratory (SML) vs. nonsulfidic male laboratory (NML)	978	377	120	57	33	17	12	3	6	4	1	9
Sulfidic female laboratory (SFL) vs. nonsulfidic female laboratory (NFL)	907	329	110	44	18	12	6	6	3	2	2	12

The genomic features associated with DMRs for each comparison were similar. The primary cytosine nucleotide followed by a guanine nucleotide (CpG) density was 1 to 6 CpG per 100 bp (Fig. 2). The average size of a DMR was around 1 kilobase (kb), with a range from 1 to 6 kb (Fig. 2). Therefore, the CpG density was low in 1 to 2 kb regions, similar to CpG deserts previously observed (47).

Fig. 2.

DMR genomic features for DMRs with a P value < 10−7. (A) The number of DMRs at different CpG densities for each of the pairwise comparisons. (B) The DMR lengths in kilobase (kb) for each of the pairwise comparisons. Sample names are: sulfidic female wild (SFW), nonsulfidic female wild (NFW), sulfidic male wild (SMW), nonsulfidic male wild (NMW), sulfidic female laboratory (SFL), nonsulfidic female laboratory (NFL), sulfidic male laboratory (SML), and nonsulfidic male laboratory (NML).

A principal component analysis (PCA) was performed to identify patterns in genome-wide methylation across all sample groups (Fig. 3). The reads per kilobase of transcript, per million mapped (RPKM) read depths at all genomic sites were used as the basis for the PCA. For samples collected in the wild, the sulfidic males and females clustered and were distinct from the nonsulfidic males and females (Fig. 3). Therefore, the sulfidic habitat influences the PCA clustering for both sexes. For the laboratory samples from fish that were propagated in a nonsulfidic environment for multiple generations, the PCA clustering again showed that the individuals derived from the sulfidic population clustered distinctly from the nonsulfidic one (Fig. 3). Hence, the samples were separated by habitat type in the PCA analysis for both the wild and laboratory populations. The PCA analysis with all samples from all comparisons further demonstrated that the wild samples clustered separately from the laboratory samples (Fig. 3). Therefore, both the population of origin and the rearing environment (wild versus laboratory) impact methylation patterns.

Fig. 3.

PCA. (A) Wild samples. (B) Laboratory samples. (C) All samples. Sample labels are: sulfidic female wild (SFW), nonsulfidic female wild (NFW), sulfidic male wild (SMW), nonsulfidic male wild (NMW), sulfidic female laboratory (SFL), nonsulfidic female laboratory (NFL), sulfidic male laboratory (SML), and nonsulfidic male laboratory (NML).

The majority of the DMRs were unique to a specific dataset at P < 10−7, with the wild or laboratory male and female comparisons having the highest level of overlap (Fig. 4). Interestingly, 94 DMRs overlapped between all four comparisons (), and nearly 20% of DMRs were shared between laboratory and field conditions at P < 10−7. We also performed an extended DMR overlap that compared the DMRs with P < 10−7 for each comparison with the genomic windows in a second comparison at P < 0.05 (Table 3). DMRs in one analysis were considered present in the second analysis if any overlapping genomic window in the second analysis had a P value of 0.05 or less. A comparison of the male and female data sets within the wild or laboratory comparisons demonstrated greater than 90% overlap. Therefore, as indicated by the PCA (Fig. 3) and Venn diagram (Fig. 4), the overlaps are very high between the sexes for each laboratory or wild samples. Interestingly, the extended overlaps for the comparisons between laboratory and wild samples for both males and females were greater than 80% (Table 3). Therefore, even after two generations of maintaining fish under nonsulfidic laboratory conditions, the majority of DMRs documented in the laboratory overlapped with those documented in the wild with the extended overlap comparison. Although this extended overlap with multiple P values is not standard, it does reveal a high degree of overlap not observed in the high statistical threshold comparison.

Fig. 4.

DMR overlap for DMRs with a P value threshold of P < 10−7. Sample groups are: sulfidic female wild (SFW), nonsulfidic female wild (NFW), sulfidic male wild (SMW), nonsulfidic male wild (NMW), sulfidic female laboratory (SFL), nonsulfidic female laboratory (NFL), sulfidic male laboratory (SML), and nonsulfidic male laboratory (NML).

Table 3.

Extended DMR Overlap

	SFW versus NFW	SMW versus NMW	SFL versus NFL	SML versus NML
SFW versus NFW		1,372 (94%)	1,207 (83%)	1,205 (82%)
SMW versus NMW	984 (94%)		871 (83%)	860 (82%)
SFL versus NFL	1,192 (82%)	1,155 (80%)		1,324 (91%)
SML versus NML	1,356 (84%)	1,291 (80%)	1,496 (92%)

The P < 10−7 left axis and P < 0.05 top axis with number of DMR overlap and percent of total for the specific comparison. Sample groups are: sulfidic female wild (SFW), nonsulfidic female wild (NFW), sulfidic male wild (SMW), nonsulfidic male wild (NMW), sulfidic female laboratory (SFL), nonsulfidic female laboratory (NFL), sulfidic male laboratory (SML), and nonsulfidic male laboratory (NML).

The DMR gene associations were investigated using DMRs identified within 10 kb of a gene so the promoter regions could be considered (). The signaling, transcription, cytoskeleton, metabolism, and receptor functional categories were the most predominant DMR-associated genes in all comparisons (Fig. 5). There were consistent similarities in the DMR-associated gene numbers for all the categories. Analysis of the DMR-associated gene functions using a Pathway Studio program analysis identified a large number of genes related to sulfur toxicity, either directly or indirectly (Fig. 6). A variety of related protein ligands, transcription factors, receptors, kinases, and phosphatases were present for the comparison data sets. Another relevant DMR-associated gene functional category involved sulfur metabolism (Fig. 6). A number of different components of sulfur metabolism processes were listed, and all had connections to DMR-associated genes. The most predominant were chondroitin sulfate–associated genes that are a critical component of the extracellular matrix and proteoglycans. The final set of DMR-associated genes analyzed was the 94 overlapping DMRs between all comparisons (Fig. 4) with the DMR list and DMR-associated genes (). These DMR-associated genes and correlated gene processes are presented in Fig. 7. Although this does not demonstrate direct links with sulfidic environments or toxicity, the pathways and gene families influenced are anticipated to be indirectly linked to the environmental exposure and potential adaptation.

Fig. 5.

DMR-associated gene categories and the number of DMRs for each comparison presented for each category. Sample groups are: sulfidic female wild (SFW), nonsulfidic female wild (NFW), sulfidic male wild (SMW), nonsulfidic male wild (NMW), sulfidic female laboratory (SFL), nonsulfidic female laboratory (NFL), sulfidic male laboratory (SML), and nonsulfidic male laboratory (NML).

Fig. 6.

DMR gene associations from Pathway Studio. (A) DMR gene associations with sulfur toxicity. (B) DMR gene associations with sulfur metabolism. Gene symbols are presented and the functional categories listed.

Fig. 7.

Overlapping 94 DMR P < 10−7 gene associations with cellular processes and pathologies from Pathway Studio. Gene symbols are presented and the functional categories listed.

Discussion

There is evidence that epigenetic modifications change in response to exposure to environmental toxicants (4, 7, 11) and that modified placement of epigenetic modifications are associated with phenotypic variation and disease (3). The analysis of methylation changes in populations adapted to a physiological stressor allows for the assessment of short-term responses to environmental stress as well as potential forecasting of long-term responses. Here, we used H2S-sensitive and H2S-tolerant fish, P. mexicana, living in naturally nonsulfidic and sulfidic environments, respectively, to elucidate the epigenetic changes in response to living in H2S-rich springs.

The current study investigates the epigenetic inheritance stability of putatively adaptive epigenetic alterations following the removal of the abiotic selective pressure. Although the P. mexicana genome (48) is not a chromosome-level assembly, sufficient information was available to assess the environmental influences on the DMRs, associated genomic features, and potential gene associations. The lack of a well-characterized genome will underestimate the DMR-gene associations, and future improved genomic information will likely improve and expand on the gene association information. The DMR size and locations in regions of low-density CpG (termed CpG deserts) found in this study are similar to environmentally responsive DMRs previously identified in other species (47, 49). Low-density CpG regions generally constitute over 90% of the genome for most species, while the high-density CpG islands generally constitute less than 5% of the genome (49). As shown in most species, from plants to mammals, CpG DNA methylation is predominant, with non-CpG methylation being negligible with no known functional significance currently identified (50). Although the MeDIP procedure used is biased to low-density CpG regions, it is the most efficient to examine the majority of the genome compared to other procedures, such as bisulfite sequencing, that can be biased to higher-density CpG analyses (47, 51). High-density CpG regions are predominantly restricted to 50% of promoters in mammals, which is less than 5% of the genome. While these regions are not detected efficiently with MeDIP sequencing, greater than 90% of the genome is detected with MeDIP sequencing (52).

The DMR-associated genes were investigated to provide insights into the potential impacts of the epigenetic alterations on the genome and physiology. Although a gene association is identified, this does not confirm a functional causal link of the epigenetic alteration and gene regulation. This is a limitation that will require future studies at the gene expression level to confirm. Due to the sulfidic-versus-nonsulfidic focus of the study, the initial analysis of the DMR-associated genes involved sulfur toxicity and metabolism. A large number of toxicity-associated genes that were either directly or indirectly correlated with sulfur that were associated with the DMRs identified including genes such as MTOR or TPMT have been shown to have direct impacts on sulfur toxicity (53–56). Sulfur metabolism genes were also associated with the DMRs identified. The chondroitin sulfate processes, which have been associated with sulfur exposures (57, 58), had the most represented associated genes. A number of cellular processes (e.g., apoptosis and angiogenesis) were also associated with the DMRs. Although these DMR-associated genes have links with sulfur toxicity and metabolism, further molecular work is needed to correlate these molecular processes and genes to adaptation of P. mexicana to H2S and the observed epigenetic inheritance.

An overlap of the DMRs between the comparisons and data sets demonstrated that the majority of DMRs were distinct for the comparison with DMRs at P < 10−7. The primary principal component in the PCA separated laboratory and wild samples, suggesting that there are largescale changes in DMRs when populations are brought into the laboratory for multiple generations. An approximately 20% overlap was observed between the male and female wild or laboratory DMR comparisons. To expand this analysis, an extended overlap was performed that compared the P < 10−7 DMR with the P < 0.05 DMRs for the other data sets. Although the extended overlap with a high and low statistical threshold is not standard, this allows for the identification of overlapping DMRs that may be marginally significant at any given P value threshold. The comparison showed a greater than 90% overlap between the male or female DMR comparisons for the wild or laboratory DMRs. Greater than 80% overlap was observed between the wild and laboratory comparisons; this level of overlap demonstrates that the majority of the sulfidic versus nonsulfidic DMRs identified were in common between the F0 generation wild populations and the F2 generation laboratory populations. Therefore, the majority of the methylation changes identified comparing wild individuals from sulfidic and nonsulfidic environments were found to be propagated for two generations in the laboratory in the absence of the sulfidic environmental exposure. These observations demonstrate an epigenetic inheritance of DMRs in the absence of the continued sulfidic environment. This is an observation of methylation inherited through the germline when an environmental stressor (i.e., H2S) is removed from a wild population and suggests that future studies on transgenerational inheritance beyond the F2 generation would be useful. This is consistent with previously described environmentally induced epigenetic transgenerational inheritance phenomena (4, 6).

Environmentally induced epigenetic inheritance of phenotypic variation and disease was first described in rats exposed to an agricultural fungicide vinclozolin (11) but now has been shown to occur in all species investigated from plants to humans (4, 6). Previous studies have primarily used the transient exposure of an individual to induce epigenetic alterations in the germline (sperm or egg) to promote the epigenetic transgenerational inheritance phenomenon (4, 6). The transmission of epigenetics between generations requires the involvement of the germline, as they are the only cells that can transmit molecular information between generations (4, 6). Direct exposure to an environmental factor can promote epigenetic alterations, and this can occur at any time during development. For example, the exposure of a gestating female mammal exposes the F0 generation mother, the F1 generation fetus, and germline that will generate the F2 generation individual (4, 59). P. mexicana is a livebearing fish that requires the same consideration. Therefore, the F3 generation is the first generation that has no direct exposure toxicity (59). For a nonlivebearing fish model, exposure of the F0 generation male or female will directly impact the F0 generation and the sperm or egg that will form the F1 generation; therefore, the F2 generation is the first generation (i.e., transgenerational) not having direct exposure toxicity (59). An example in fish involves the exposure of zebrafish to mercury that promoted neurodevelopmental abnormalities in zebrafish in the F2 generation offspring, which was the first epigenetic transgenerational inheritance observation in fish (60). However, livebearing fish are more like mammals in that the exposure of a gestating female will expose the germline that will generate the F2 generation. The current study observations suggest that an important future study in these livebearing fish is the investigation of DMRs in the F3 generation, the first with no direct exposure to H2S, which would be considered epigenetic transgenerational inheritance. Although the current study did not analyze the F3 generation in the nonsulfidic environment, propagation of the F2 generation in the laboratory environment demonstrates survival of the fish to the F3 generation and beyond. Future studies will need to compare the F1 generation of laboratory fish to determine whether there is a steady decrease in DMR shared with each subsequent generation or whether the changes are due to the initial transition to an aquarium environment and are stable and maintained at similar levels thereafter. It is also unknown whether fish from sulfidic habitats that have been raised in a nonsulfidic laboratory environment for multiple generations may have a higher survival than fish with nonsulfidic ancestry when transmitted back into a sulfidic environment. The question is whether an ancestral sulfidic environment adapted epigenetics is maintained, which may facilitate reintroduction into the sulfidic environment. This will be an important future experiment to consider for understanding how epimutations altered in an aquarium environment play a role in phenotypic variation and survival.

Although the ability of environmental stressors to alter the epigenome of the individual exposed and impact the physiology of the individual is critical, the ability to alter the germline epigenome and propagate this transgenerationally to subsequent generations allows for an adaptive mechanism that can impact evolution (61). The phenotypic plasticity and molecular alterations in the epigenome and transcriptome to survive the sulfidic environment are essential to allow an organism to respond to its environment. The epimutations documented here may contribute to heritable variation in gene expression, which we have documented in detail in the field, the laboratory, and across multiple species (40, 41, 44, 62). The sulfidic and nonsulfidic populations of P. mexicana studied here are genetically differentiated, with a divergence time of ∼10,000 y ago and an average empirical FST of 0.11 (63), but considered the same species (37–39), and studying how these epigenetic changes play into adaptive changes will be a fruitful avenue of research. The population studied here is one of four pairs of sulfidic and nonsulfidic P. mexicana populations in tributaries of the Río Grijalva, and the replicated nature of the independently derived sulfidic populations provides an opportunity to test whether the findings in this drainage are consistent across multiple drainages. Future research will need to explore how genetic variation interacts with epimutations to shape phenotypic variation in this system. This project is a first step to evaluate epigenetic effects of toxicant exposure and sets the stage for future studies examining how long-term environmental exposure may shape the epigenome.

Materials and Methods

Samples.

Male and female P. mexicana were collected from H2S-rich habitats and nonsulfidic habitats in southern Mexico (Fig. 1) for the preparation of purified RBCs. Wild-caught samples were obtained from southern Mexico in proximate H2S-rich (sulfidic, n = 4 female, n = 4 male) (PSO, Latitude: 17.43843, Longitude: −92.77476) and nonsulfidic springs (n = 4 female, n = 3 male) (Bonita, Latitude: 17.42685, Longitude: −92.75213) in the Tacotalpa Drainage in southern Mexico. Another set of P. mexicana males and females from sulfidic or nonsulfidic habitats were collected and transported into a breeding population at Kansas State University and propagated for two generations (F2) in the nonsulfidic laboratory environment (Fig. 1) (sulfidic, n = 4 female, n = 4 male) and (nonsulfidic, n = 4 female, n = 4 male). All laboratory animals were reared in nonsulfidic water, and n = 4 individuals were sampled per population. Adult fish for meDIP sequencing were collected in the wild in May 2015. For laboratory experiments, pregnant females were collected in May 2013, and adult F2 offspring were harvested in August 2015 for blood extractions. Of note, fish from sulfidic habitats are smaller than fish from the nonsulfidic habitats, and males are smaller than females.

DNA Extraction.

Fish were euthanized by cervical transection, and blood was immediately extracted by using microhematocrit capillary tubes (Fisher). Tubes were then spun down at 12,000 rpms for 3 min using a ZipCombo Centrifuge (LW Scientific) to separate blood cells from serum. After centrifugation, capillaries were broken to retain only RBCs. The portion of the capillary containing the RBCs was placed in a 1.5-mL Eppendorf tube for DNA extraction. DNA was extracted using the Qiagen DNeasy Blood & Tissue kit according to the manufacturer’s instructions.

MeDIP.

Extracted genomic DNA was used to isolate methylated DNA with a methyl-cytosine antibody precipitation procedure (MeDIP). The protocol is described in detail in ref. 64. A total of 5 μg total genomic DNA was sonicated using a Covaris sonicator. Sonicated DNA was diluted with TE (tris-ethylenediaminetetraacetic acid [EDTA]) buffer to 400 μL, heat denatured for 10 min at 95 °C, and cooled on ice for 10 min to create single-stranded DNA fragments. A total of 100 μL of 5× IP (immunoprecipitation) buffer and 5 μg of antibody (monoclonal mouse anti 5-methyl cytidine; Diagenode C15200006) were added to the fragmented single-stranded DNA. The mixture was incubated on a rotator overnight at 4 °C. Prewashed magnetic beads (50 μL, Dynabeads M-280 Sheep anti-Mouse IgG; Life Technologies 11201D) were added to the 500 μL of DNA–antibody mixture. The mixture was incubated for 2 h on a rotator at 4 °C, after which the samples were washed three times with 1× IP buffer. Washed samples were resuspended in 250 μL digestion buffer (5mM Tris pH 8, 10.mM EDTA, and 0.5% SDS [sodium dodecyl sulfate]) with 3.5 μL Proteinase K (20mg/mL) and incubated at 55° for three hours. A Phenol-Chloroform-Isoamyalcohol extraction was used to clean-up the DNA, and the DNA was resuspended in 20 μL of H2O. DNA concentration was measured in Qubit (Life Technologies) with the single-stranded (ss) DNA kit (Molecular Probes Q10212).

MeDIP Sequencing.

For library preparation, we used the NEBNext Ultra RNA Library Prep Kit for Illumina starting at step 1.4 of the manufacturer’s protocol to generate double-stranded DNA. After this step, the manufacturer’s protocol was followed. Each individual received a unique barcode. Enriched methylated DNA libraries were sequenced at the Washington State University Genomics Core in Spokane, WA, using an Illumina HiSeq 2500 with paired-end 50-bp reads.

Identifying and Analyzing Differentially Methylated Regions.

DMR identification and annotation followed previously published aapproaches (65, 66). The FastQC program (67) was used to assess data quality. Reads were trimmed to remove adapters and low-quality bases using Trimmomatic (68). The reads for each MeDIP sample were mapped to the P. mexciana (48) genome using Bowtie2 (69) with default parameter options. The mapped read files were then converted to sorted BAM (Binary Sequence Alignment/Map) files using SAMtools (70).

Differential coverage between sulfidic and nonsulfidic populations was calculated using the MEDIPS R package (71). P value from edgeR (72) was used to determine the significance of the difference between the two groups for each 100-bp genomic window. Windows with an edgeR P value less than a specified threshold (P < 10−7) were considered the initial start of the DMR. DMR edges were extended until no genomic window with a P value less than 0.1 remained within 1,000 bp of the DMR. The extended DMR overlap compared the DMRs, with at least one 100-bp window with P < 10−7 from one comparison, with the genomic windows (100-bp regions) in a second comparison. Windows that had a P value <0.05 in the second comparison were considered overlapping. The Ensembl database (73), accessed with the biomaRt R package (74), was used to annotate DMRs. Genes that were overlapping a DMR, including 10 kb on either side of the DMR, were input into a Kyoto Encyclopedia of Genes and Genomes pathway search (75, 76) to identify relevant associated pathways. The DMR-associated genes were sorted into functional groups using information provided by the DAVID (Database for Annotation, Visualization and Integrated Discovery) (77) and Panther (78) databases incorporated into an internal curated database (https://www.skinner.wsu.edu/). DMR-associated gene correlations present in published literature were further analyzed using Pathway Studio software (version 12.2.1.2: Database of functional relationships and pathways of mammalian proteins; Elsevier). R code computational tools are available at GitHub (https://github.com/skinnerlab/MeDIP-seq) and https://www.skinner.wsu.edu.

Ethical Statement.

Procedures for all experiments were approved by the Institutional Animal Care and Use Committee at Kansas State University (protocol no. 3418). Field work was approved by the Mexican government (fieldwork permit DGOPA.00093.120110.-0018).