Abstract
Epigenetic inheritance has been proposed to contribute to adaptation and acclimation via two information channels: (i) inducible epigenetic marks that enable transgenerational plasticity and (ii) noninducible epigenetic marks resulting from random epimutations shaped by selection. We studied both postulated channels by sequencing methylomes and genomes of Baltic three-spined sticklebacks (Gasterosteus aculeatus) along a salinity cline. Wild populations differing in salinity tolerance revealed differential methylation (pop-DMS) at genes enriched for osmoregulatory processes. A two-generation experiment demonstrated that 62% of these pop-DMS were noninducible by salinity manipulation, suggesting that they are the result of either direct selection or associated genomic divergence at cis- or trans-regulatory sites. Two-thirds of the remaining inducible pop-DMS increased in similarity to patterns detected in wild populations from corresponding salinities. The level of similarity accentuated over consecutive generations, indicating a mechanism of transgenerational plasticity. While we can attribute natural DNA methylation patterns to the two information channels, their interplay with genomic variation in salinity adaptation is still unresolved.
INTRODUCTION
Recent advances in epigenetics challenge our understanding of inheritance and adaptive evolution (1–3). It has been suggested that epigenetic modifications—for example, via DNA methylation, histone modification, or small RNAs—create phenotypic diversity and ultimately contribute to rapid evolutionary adaptation (4–6). Several theoretical models posit that the heritable proportion of these molecular modifications can be classified into two distinct information channels (5, 7, 8). Selection-based epigenetic marks emerge as spontaneous epimutations that remain stable across subsequent generations, although their overall stability is three to four orders of magnitude lower compared to DNA base changes (7, 9). Similar to adaptation from DNA sequence–based variation, these epimutations may result in different phenotypes that become targets of natural selection and thereby carry information on past selection regimes without directly responding to the current environment (5, 8, 10). On the other hand, detection-based effects describe inducible epigenetic marks at defined genomic locations, which are under environmental control (7). Such transfer of parental information linked to environmental cues represents a rapid and reliable mechanism underlying transgenerational plasticity, which is hypothesized to buffer the extinction risk of populations under sudden environmental change until genetic adaptations can catch up (“genetic rescue”) (7, 11). Distinguishing between these mechanisms is important because they have very different implications for the evolution of populations. Stable epigenetic marks follow evolutionary principles of DNA sequence–based inheritance with random variation shaped by selection. In contrast, directional processes via inducible epigenetic marks can be considered a transgenerational form of plasticity that involve previously evolved regulatory mechanisms targeting specific sites on the genome. While the principal differences between these two transmission channels are clear (4, 5, 7), empirical evidence for their presence in wild vertebrate populations is lacking.
Here, we assess whether these two epigenetic information channels can be detected in nature and test whether short-term acclimation responses match patterns of DNA methylation variation of locally adapted populations. Transgenerational experiments that yield DNA methylation profiles more similar to those of locally adapted natural populations would provide evidence that DNA methylation is mechanistically involved in adaptive transgenerational plasticity.
Studying adaptation to ocean salinity is particularly suited to identification of selection- and detection-based effects because spatiotemporal patterns in ocean salinity are more stable than other variables, for instance, temperature. Since salinity change imposes strong physiological stress with well-defined cellular effects (12), natural salinity gradients offer unparalleled opportunities to use local patterns of epigenetic variation as background against which direction and magnitude of results from experimental salinity manipulations can be tested. One suitable ecosystem to follow such a space-for-time approach is the Baltic Sea, which is a semi-enclosed marginal sea that has been dubbed a “time machine” to evaluate the predicted perturbations associated with global change (13).
Taking advantage of the Baltic Sea salinity gradient, we sequenced the methylomes [reduced representation bisulfite sequencing (RRBS)] and whole genomes of three-spined sticklebacks (Gasterosteus aculeatus) from three populations that are locally adapted to different salinities [6, 20, and 33 practical salinity units (PSU)] (14) in and outside the Baltic Sea. Specifically, we focus on the patterns of (epi)genomic variation, while transgenerational phenotypic effects have been described previously for the exact same populations (15).
Baltic stickleback populations are genetically differentiated [genome-wide average pairwise FST = 0.028 (14)] and show patterns consistent with local adaptation to salinity regimes in controlled common garden experiments (15, 16). Moreover, previous studies have revealed transgenerational plasticity in response to variation in temperature (17) and changes in DNA methylation levels at osmoregulatory genes in response to within-generational salinity manipulation (18, 19). However, it remains unclear whether DNA methylation mediates transgenerational plasticity, a possible mechanism enabling adaptive phenotypes to rapidly emerge in the face of environmental change. In this study, we consider transgenerational effects to be adaptive if the preacclimation of the parents enhances the fitness of the offspring, sometimes referred to as intergenerational effects. To address this question, we complemented our field survey with a two-generation salinity acclimation experiment using the mid-salinity population (20 PSU). This experiment enabled us to quantify the proportion of noninducible (stable, potentially selection-based) and inducible (potentially detection-based) DNA methylation within and across generations (Fig. 1), acknowledging that we tested methylation mark stability only with respect to experimental salinity manipulation. We focused on the methylation of cytosines at cytosine-phosphate-guanine dinucleotides (CpG sites), the most common methylation motif in vertebrates (20), with partial inheritance potentially involved in adaptive evolution (11).
We tested three nonexclusive hypotheses: (i) Stickleback populations from different salinities (6, 20, and 33 PSU) show differentially methylated CpG sites (hereafter referred to as pop-DMS). (ii) Such pop-DMS include both types of methylation sites: experimentally stable sites (potentially selection based) and experimentally inducible sites (potentially detection based). (iii) Upon transgenerational salinity acclimation, inducible DNA methylations become more similar to the patterns of natural populations at corresponding salinities. When associated with beneficial phenotypic effects and increased relative fitness, the latter would be evidence for a mechanism underlying adaptive transgenerational plasticity (overview in Fig. 2).
RESULTS
Identifying differentially methylated CpG sites between stickleback populations along a natural salinity cline
pop-DMS were determined via RRBS in 46 wild-caught sticklebacks from three different sites that varied in average salinity [Sylt (SYL), 33 PSU; Kiel (KIE), 20 PSU; Nynäshamn (NYN), 6 PSU; Fig. 1]. After quality and coverage filtering, we obtained 525,985 CpG sites present in all groups (q < 0.0125; methylation difference, ≥15%), corresponding to ~4% of all CpG sites in the stickleback genome. Among pairs of wild-caught populations, we detected 1470 (comparison of 20 versus 6 PSU) and 1158 (20 versus 33 PSU) pop-DMS. The distribution of these sites was random with regard to the genomic features (promoter, exon, intron, and intergenic; 20 versus 6 PSU: X23 = 3.36, P = 0.340; 20 versus 33 PSU: X23 = 1.61, P = 0.656; table S1) and chromosomal regions (fig. S1A). Among these pop-DMS, 1098 (20 versus 6 PSU) and 871 (20 versus 33 PSU) were located close to [<10 kb from transcription start sites (TSS)] or within genes and thereby associated with 655 and 510 genes, respectively. Many of these genes are involved in fundamental biological processes such as DNA repair and strand renaturation, as well as chromosome condensation and separation (fig. S2). Of particular relevance is the enrichment in genes associated with osmoregulatory processes such as ion transport and channel activity, renal water homeostasis and absorption, and urine volume regulation (Fig. 3). Genes associated with ≥10 pop-DMS are listed in Table 1 [for all genes, see table S2 (A and B)].
Characterizing stable and inducible DNA methylation in a two-generation salinity acclimation experiment
To assess the proportion of inducible DNA methylation, we conducted a two-generation salinity acclimation experiment with laboratory-bred sticklebacks from the mid-salinity population that was subjected to either increased or decreased salinity (Fig. 1). We considered pop-DMS to be noninducible (hereafter referred to as “stable”) when both the within-generational and the transgenerational acclimation groups were not differentially methylated compared to the control group (q ≥ 0.0125). On the other hand, if a pop-DMS was differentially methylated between at least one of the acclimation groups (within- and transgenerational) compared to the control group (q < 0.0125; methylation difference, ≥15%), then this site was considered inducible. Pop-DMS with a significant q value not exceeding the threshold of differential DNA methylation were treated as a separate category (hereafter referred to as inconclusive). After two generations of salinity acclimation, we found that most of the pop-DMS remained stable, regardless of the direction of salinity change (926 pop-DMS, 63% at decreased salinity; 694 pop-DMS, 60% at increased salinity). A smaller number of pop-DMS (13%) were inducible, as they showed a significant change in CpG methylation upon experimental salinity decrease (198 pop-DMS) or increase (148 pop-DMS). An additional 24 and 27% (346 and 316 pop-DMS, respectively) were inconclusive. The number of inducible pop-DMS (13%) derived from comparisons between natural populations was much higher than expected from a random subset of CpG sites across the genome (<1%; 1000 replicates; salinity decrease: X22 = 1090.7, P < 0.001; salinity increase: X22 = 967.7, P < 0.001). This means that pop-DMS are enriched for sites that plastically respond to salinity change, which is expected for populations from different salinities.
Stable and inducible pop-DMS are associated with different functional gene categories
Gene functions associated with stable pop-DMS (452 and 329 under salinity decrease and increase, respectively) were enriched not only for a number of fundamental biological processes such as DNA repair and chromosome separation (fig. S2) but also for osmoregulatory functions (e.g., ion channel activity; Fig. 3). Furthermore, under increased salinity, many metabolic processes were found among the stable pop-DMS (Fig. 3). Inducible pop-DMS were associated with genes (100 and 82 under salinity decrease and increase, respectively) that were primarily enriched for other osmoregulatory functions regulating, for example, ion transmembrane transport (Fig. 3 and fig. S2). Therefore, stable and inducible pop-DMS affect not only different genes but also different gene ontologies with little overlap (Fig. 3 and fig. S2).
Assessing the role of inducible DNA methylation in nature
We investigated whether multiple pop-DMS associated with the same gene showed a correlated response to experimental salinity acclimation, which would require that they are nonrandomly distributed among the three categories stable, “inducible,” and “inconclusive.” Accordingly, we found a correlated response for pop-DMS at 13 of 20 genes (genes with more than 10 pop-DMS; Fisher’s exact test, P < 0.05; Table 1), which suggests that inducible pop-DMS are predefined and directed.
We then tested whether inducible pop-DMS in the experimental fish became more similar to methylation levels found in natural populations. Of the 198 (decreased salinity) and 148 (increased salinity) inducible pop-DMS, 130 (66%) and 101 (68%), respectively, became more similar to methylation levels of wild population to the corresponding salinity (hereafter referred to as “expected” direction). Conversely, at 68 (34%; decreased salinity) and 47 (32%; increased salinity) inducible pop-DMS, experimental fish showed methylation changes in the opposite direction, reducing the similarity to methylation levels observed in the natural populations (hereafter referred to as “opposite” direction).
Why, in a proportion of inducible methylation marks, the similarity between experimental and natural methylation levels was reduced was puzzling. One explanation could be a high level of genomic differentiation between the populations at these sites since genomic variation can have a strong cis-regulatory impact on epigenomic variation and may alter direction and function of methylation marks together (21). Thus, we hypothesized that opposite inducible pop-DMS are more often occurring in regions with higher genomic (DNA sequence–based) differentiation, while we anticipated the reverse at expected inducible pop-DMS. Accordingly, we resequenced whole genomes of the same wild-caught individuals that we used for RRBS and calculated the degree of genomic differentiation per inducible pop-DMS as mean FST value (±5-kb window) between populations. In line with our hypothesis, the populations from KIE (20 PSU) and NYN (6 PSU) were genetically more differentiated at opposite inducible pop-DMS than at expected sites (decreased salinity: δ.mean.FST = −0.014, P = 0.002; Fig. 4, A and C). A similar, yet not significant, trend was found between the populations from KIE (20 PSU) and SYL (33 PSU) (increased salinity: δ.mean.FST = −0.005, P = 0.153; Fig. 4, B and D). An alternative explanation is that not only salinity but also, rather, a combination of environmental cues (i.e., temperature, predation, and food) resulted in the methylation patterns found in the SYL population, which we did not include in our experiment. To understand whether selection has shaped the differences between increased and decreased salinity exposure, we tracked survival rates from fertilized eggs to the 3-month-old offspring and compared them between treatment groups. Mortality differed significantly between the treatment groups [generalized linear mixed model (GLMM), X24 = 66.159, P < 0.001; Fig. 5A and table S3A] with increased mortality under increased salinity, while mortality under decreased salinity was generally low and did not differ from the control group (Fig 5A and table S3A). Hence, while we cannot entirely disregard the effect of selection for increased experimental salinity, the patterns observed at pop-DMS upon reduced salinity are likely the sole result of tolerance mechanisms for salinity change.
Comparing within- and transgenerational acclimation effects on inducible DNA methylation
To test for adaptive transgenerational plasticity, we evaluated whether salinity acclimation over two, instead of only one, consecutive generations enhances the similarity of inducible pop-DMS with patterns found among wild populations at corresponding salinities. To do so, we calculated the percentage match (δ.meth.diff; Fig. 4, E to H) between the experimental groups and the anticipated adaptive methylation levels of wild populations. In line with our hypothesis, we found that transgenerational compared to within-generational salinity manipulation increased the δ.meth.diff (for expected inducible methylation, decreased salinity: F1,256 = 30.42, P < 0.001; increased salinity: F1,198 = 10.39, P = 0.001; Fig. 4, E and F). Under decreased experimental salinity, we found an interaction of “methylation direction” (hyper- or hypomethylation) and “acclimation” (within- and transgenerational) affecting the δ.meth.diff [analysis of variance (ANOVA), δ.meth.diff ~ methylation direction × acclimation, F1,256 = 7.69, P = 0.006; Fig. 4E]. Specifically, transgenerational acclimation increased the similarity of hypomethylated sites to methylation levels found in natural populations, while hypermethylated sites showed similar values within and across generations (Fig. 4E). While, for expected inducible sites, this effect was only present under decreased salinity, at opposite inducible sites, transgenerational acclimation to decreased and increased salinity elevated the δ.meth.diff at hypomethylated sites (ANOVA, δ.meth.diff ~ methylation direction × acclimation, decreased salinity: F1,132 = 19.89, P < 0.001; increased salinity: F1,90 = 9.85, P = 0.002; Fig. 4, G and H).
To infer the effect of DNA methylation differences on offspring, we compared fitness proxies among control, within-generational, and transgenerational acclimation groups (22). Specifically, we assessed the total weight, standard length (SDL), and the hepatosomatic index (HSI) as a proxy for energy reserves in the form of liver glycogen storage. SDL (GLMM, X24 = 9.965, P = 0.041; Fig. 5B and table S3B) and total weight (GLMM, X24 = 11.518, P = 0.021; Fig. 5D and table S3D) differed between treatment groups. Highly significant differences were detected for the HSI (GLMM, X24 = 22.688, P < 0.001; Fig. 5C and table S3C), with elevated HSI observed under decreased salinity compared to fish from the control group. This supports previous findings, showing that osmoregulation at 6 PSU is energetically less demanding than that at higher salinities (15, 16). Under increased salinity, HSI was lower compared to that of fish exposed to decreased salinity in the within-generational acclimation group, while a transgenerational acclimation to increased salinity partially removed this difference. Although not significant, we observed a trend toward higher mean HSI in the transgenerational acclimation group compared to the within-generational acclimation group at the same salinity (Fig. 5C and table S3C).
DISCUSSION
This study investigated whether two postulated channels of epigenetic inheritance (selection based and detection based) can be identified in natural populations, focusing on salinity adaptation among populations of three-spined sticklebacks. Consistent with expectations for selection-based DNA methylation sites (7), we identified pop-DMS between populations that were both enriched for osmoregulatory functions and stable with respect to two generations of experimental salinity manipulation. Phenotypic variation originating from selection-based DNA methylation sites that are expected to have high epimutation rates [~10−4 for Arabidopsis thaliana (9)] could allow populations to explore the fitness landscape faster than under DNA sequence–based genetic variation alone [mutation rate, ~10−8 (5, 23)]. Whether and at which rate these randomly emerging epimutations, as predicted for selection-based DNA methylation (5), occur in vertebrates remains unresolved. Notwithstanding, the observed enrichment of osmoregulatory gene functions for stable methylation sites (Fig. 3) suggests that they were subject to divergent natural selection, possible in interaction with DNA sequence–based variation. Furthermore, since local adaptation is 10 times more likely to involve changes in gene expression than in amino acid sequence (24), it is conceivable that differential DNA methylation and, consequently, regulation of osmoregulatory genes may contribute to local salinity adaptation. In sticklebacks, for instance, immunological adaptation has been shown to be mediated by gene expression (25). One of the top candidate genes differentially methylated between populations from 20 and 6 PSU was eda (ectodysplasin A), a well-described gene involved in lateral plate formation (26). Salinity and calcium are significant drivers of plate morphology (27) in proposed conjunction with predation (28). Our findings suggest that repeated and parallel selection for the low plated eda allele in response to low saline habitats (29–31), including the Baltic Sea (14, 32), may also involve methylation-related mechanisms. Previous studies have shown that energetic cost for Baltic sticklebacks increases with increasing difference between treatment and isosmotic salinity conditions [~11 PSU (33)] (15, 16). In line with these findings, we observed many metabolic processes associated with stable pop-DMS under increased salinity, also reflected in the lower HSI of fish at that salinity. Together, our results on the noninducible fraction of differentially methylated genes are consistent with a role in local salinity adaptation across stickleback populations (Fig. 3; Table 1; fig. S2; and table S2, A and B). These patterns of local adaptation in DNA methylation can have a genomic basis in the form of cis- and trans-acting genomic loci (21, 34). Whether the differential methylation patterns represent an independent mechanism for local adaptation or are rather a consequence of DNA sequence–based genetic differentiation needs further study. Since our experiment only manipulated salinity while keeping all other factors constant, it is possible that some pop-DMS that were stable under salinity change could be inducible by other changing parameters.
With respect to the second postulated information channel, detection-based epigenetic inheritance (7), we identified more experimentally inducible pop-DMS than expected by chance. Multiple DMS associated with the same gene showed synchronized responses (Table 1). Furthermore, inducible pop-DMS were associated with different osmoregulatory genes compared to stable pop-DMS. Thus, inducible sites reflect a salinity-mediated plastic response, allowing individuals to regulate their ion balance relative to the seawater medium instantaneously without requiring any further genetic adaptation. More than two-thirds of these inducible pop-DMS became more similar to methylation patterns found in wild population. The similarity of these pop-DMS methylation levels between naturally adapted and experimentally acclimated population increased across generations. Considering the corresponding beneficial phenotypic effects, this strongly suggests that adaptive transgenerational plasticity plays a role in salinity acclimation. Since we used a split-clutch design for the breeding experiment, we can assume that these groups have similar genomic backgrounds. Furthermore, as mortality levels at low salinity remained low and did not differ between treatment groups, we can rule out any effect of selection altering the genotype composition in the groups at decreased salinity.
The induction of methylation sites has been discussed as a potential buffer for environmental changes (11, 17, 35). We found that the potential for adaptive transgenerational effects, specifically the ability to establish the anticipated adaptive methylation pattern found in the wild, differed among methylation directions (hypo- and hypermethylated sites; Fig. 4, E, G, and H), with a higher potential for transgenerational plasticity at hypomethylated sites. In line with our finding, the spontaneous addition of a methyl group to a cytosine is 2.5 times more likely than the removal (23). Methylation reprogramming that includes extensive methylation removal and de novo methylation during gamete formation and zygote development could thus serve as mechanisms to demethylate CpG sites in the transgenerational acclimation group (36, 37).
The genetic background is considered to be an important source for epigenomic variation via cis- and trans-regulatory mechanisms (21, 38, 39). Thus, we characterized the genomic region surrounding each inducible pop-DMS and quantified the level of population differentiation (FST). This analysis revealed a negative correlation between population genetic differentiation and the propensity of the experimental population to approach the methylation level of the low salinity population (NYN) under salinity decrease (Fig. 4, A and C). Here, experimentally induced DNA methylation becomes more similar to the methylation in natural populations only in genomic regions with low genetic differentiation. On the other hand, when experimentally induced methylation differences to the low salinity population increase (Fig. 4, A and C), this occurs in a more divergent genomic background, suggesting that the genome has undergone selection leading to DNA-based local adaptation, rendering epigenetic modifications less relevant (5). Under increased salinity, a relationship between genomic differentiation (as FST) and methylation direction was inconclusive, suggesting that a combination of environmental cues shaped DNA methylation levels among wild populations at these sites. Overall, these findings emphasize the importance of the genomic background for interpreting DNA methylation patterns.
Together, our study provides the first empirical evidence that stable and inducible DNA methylation in wild animal populations follows predictions from evolutionary theory of selection- and detection-based epigenetic information channels (Fig. 2) (5, 7). While the selection-based information channel assumes random variation from epimutation that is subsequently shaped by selection or drift, the detection-based information channel allows a directional response in the form of transgenerational plasticity. Because the evolutionary implications of these two channels of inheritance are very different, future transgenerational or epigenetic studies should distinguish among both fundamentally different processes. Whether epigenetic marks, such as differentially methylated sites studied here, can permanently be attributed to one of the two categories or rather represent a continuum of stability levels and directionality will need further experimental testing over multiple generations.
MATERIALS AND METHODS
Animal welfare
All catches were performed under legal authorization issued by the German Ministry of Energy Transition, Agriculture, Environment, Nature and Digitalization in Schleswig-Holstein (MELUR: V242-7224.121-19), the Danish Ministry of Food, Agriculture and Fisheries of Denmark (case no: 14-7410-000227), the Estonian Ministry of the Environment (Keskkonnaministeerium – eripüügiluba nr 28/2014), and the Swedish Sea and Water Authority (Havs och Vattenmyndigheten). Ethical permission for the experiments required by German law was given by the MELUR: V312-7224.121-19, and the study is also in line with the Institutional Animal Care and Use Committee guidelines.
Survey and experimental design
For the field survey, we collected juvenile three-spined sticklebacks (G. aculeatus; 31.68 ± 14.25 mm) from three different salinity regimes inside and outside the Baltic Sea [SYL, Germany (55°00′58.3″N, 8°26′22.0″E), 33 PSU (n = 16); KIE, Germany (54°26′11.8″N, 10°10′20.2″E), 20 PSU (n = 16); NYN, Sweden (58°52′44.7″N, 17°56′06.2″E), 6 PSU (n = 16)] in September 2014. Fish were immediately euthanized using tricaine methane sulfonate solution (MS222), photographed, measured (length and total weight), and stored in RNAlater solution (24 hours at 7°C, afterward at −20°C). A cut along the ventral side ensured that the RNAlater solution diffused into all tissues. Conserved specimens were later dissected in the laboratory, and gill tissue was separated as the main osmoregulatory organ in fishes. For the acclimation experiment, we collected adult fish from KIE (20 PSU), which were crossed in our facilities at GEOMAR to obtain 10 F1 laboratory-bred families, herein referred to as “parental generation”. At 9 months after hatch, we split each family into three salinity treatment groups of 10 fish each: one at 33 PSU, one at 6 PSU, and one control group at 20 PSU. The salinity transition was performed within 10 days by 3-PSU steps every second day. Over the entire time, each group was fed ad libitum and kept in a 20-liter aquarium connected to one of three filter tanks per salinity treatment. After 5 months under treatment conditions, six pure crosses per salinity treatment group were performed in vitro, herein referred to as “offspring generation” (F2). Offspring and parental generations were kept at 18°C water temperature and a 15:9 light/dark (L/D) cycle. During the past 8 weeks before the F2 crosses, the F1 generation underwent an artificial winter to trigger reproduction (2 weeks at 12°C, 12:12 L/D; 4 weeks at 6°C, 8:16 L/D; 2 weeks at 12°C, 12:12 L/D). Upon fertilization, clutches were split and separated into different treatments (Fig. 1). At 3 months after hatch, laboratory-bred F2 sticklebacks were euthanized using MS222, photographed, and dissected, and their gill tissue was stored in RNAlater solution. The age at sampling matched the estimated age of the wild-caught juveniles (3 months). In addition to the 48 wild-caught individuals from KIE, NYN, and SYL that were used in the above field survey, we sequenced whole genomes from gill tissue of an additional three populations of sticklebacks, namely, from Falsterbo, Sweden (55°24′46.6″N, 12°55′52.3″E; 10 PSU; n = 16), Letipea (59°33′07.6”N, 26°36′29.7″E; 4 PSU; n = 16), and Barsta (62°51′47.1″N, 18°23′51.0″E; 5 PSU; n = 16).
Mortality and HSI
Mortality was monitored throughout the experiment to account for possible nonrandom effect of selection. Three months after hatch, we assessed the SDL, total weight, and liver weight of the experimental F2 generation and calculated the HSI (HSI = liver weight/total weight × 100), which is a proxy for energy reserves in the form of glycogen storage. We analyzed the effect of treatment (five treatment groups; Fig. 1) on the survival rate per family as a ratio of “alive” versus “dead” fish using glmer implemented in the R package “lme4” (40) with binomial error and “crossing” as well as “climate chamber” as random effects. The effect of treatment on HSI, SDL, and total weight was analyzed fitting three individual linear mixed-effect models using lmer in lme4 (40) with Gaussian error and crossing as well as tank nested within climate chamber as random effects. Tukey post hoc tests were run using the glht function implemented in the package multicomp (41) to identify significant differences between treatment groups.
DNA extraction
For the field survey, DNA extraction of gill tissue (n = 16 individuals per population) was performed using the DNeasy Blood and Tissue Kit (QIAGEN). Further purification of the extracted DNA was done with NucleoSpin gDNA Clean-up (Macherey-Nagel). For laboratory-bred F2 offspring of the two-generation acclimation experiment, dual extraction of whole RNA and DNA was performed from gill tissue (n = 11 to 12 individuals per treatment group; Fig. 1) stored in RNAlater solution using the AllPrep DNA/RNA Mini Kit (QIAGEN). Purity and quality of the extracted DNA were estimated using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific) and a standard agarose gel (1% agarose/tris-acetate-EDTA). DNA concentration was assessed using the Qubit 2.0 Fluorometer (Thermo Fisher Scientific). To obtain a balanced sex ratio (50:50), we determined the gender of the individuals using a sex-specific genetic polymorphism in isocitrate dehydrogenase with a modified protocol from Peichel et al. (42). For the polymerase chain reaction (PCR) (settings: once 94°C for 3 min; 30 cycles of 94°C for 30 s, 54°C for 20 s, and 72°C for 30 s; once 72°C for 5 min), 1 μl of forward and reverse primer (5 μM) was used with 4.9 μl of water, 1 μl of 10× buffer, 1 μl of deoxynucleotide triphosphate (0.5 μM), and 0.1 μl of DreamTaq (5 U/μl). The resulting PCR products were visualized with a capillary electrophoresis on the 3100 ABI sequencer and a 500 LIZ size standard. While males show a heterogametic signal with two bands [at approximately 300 and 270 base pairs (bp)], females lack the band at 270 bp.
Library preparation and sequencing (whole-genome sequencing)
For whole-genome sequencing (WGS), the TruSeq Nano DNA (Illumina) library preparation kit was used according to the manufacturer’s protocol by the Sequencing Facility of the IKMB, University of Kiel. Ultrasonication was conducted with a Covaris E220 (Covaris) to shear the input DNA (100 ng per sample and 350-bp insert size). Before the enrichment with a PCR step (8 cycles), fragmented and bead-purified DNA was ligated with adenylate at the blunt 3′ ends (end repair and A-tailing) and with indexing adapters. Fragments were cleaned with MagSi-NGSprep Plus Beads (Steinbrenner). Paired-end sequencing of the quality-controlled and multiplexed libraries was performed on the Illumina HiSeq 4000 platform (2 × 150–bp reads).
Quality assessment, data filtering, and mapping (WGS)
The command line tools of Picard version 2.7.1 (Broad Institute 2016) was used to (i) reformat the Fastq to uBAM file format and to add further values (read group, etc.) to the SAM header using FastqToSam, (ii) mark the location of adapter sequences using MarkIlluminaAdapters, and (iii) reconvert the sequences to Fastq format with SamToFastq. The stickleback genome (Broad/gasAcu1) was indexed with bwa index and used as a reference for the mapping with bwa mem (43) version 07.12-r1044. To retain the meta-information from the uBAMs, we used MergeBamAlignment. Picard was also used to identify duplicates with MarkDuplicates. Basic statistics were generated with CollectWgsMetrics, CollectInsertSizeMetrics, and AlignmentSummaryMetrics and summarized with MultiQC version 1.0.dev0 (44). A total number of 4,463,070,154 high-quality reads (mapping quality, >Q20) was mapped resulting in a mean depth of 13.84× (sd. 2.02×) and a mean insert size of 383.07 bp (sd. 9.40 bp; table S3). GATK version 3.7 HaplotypeCaller (45) was run to determine the likelihoods of the haplotypes per sample, i.e., to call single-nucleotide polymorphisms (SNPs) and indels (insertion-deletion), which were then processed with GenotypeGVCFs for a joint genotyping. SNPs were selected using hard filters for quality and extracted from the raw genotypes with a combination of the SelectVariants, VariantsToTable, and VariantFiltration commands. VCFtools (46) was used in a next step, removing SNPs with a minimum quality score below 20 and a minor allele frequency greater than or equal to 0.0049.
Library preparation and sequencing (RRBS)
The library preparation for methylation analyses followed the Smallwood and Kelsey RRBS protocol (47). A total of 100- to 250-ng purified DNA was digested with the methylation-insensitive Msp I restriction enzyme, which cuts at the “CCGG” motif and thereby enriches for CpG regions. DNA end-repair and A-tailing were conducted, and untailed CEGX spike-in controls (Cambridge Epigenetix) were added. These are DNA oligos of known sequence and with known cytosine modification, which can be used for downstream assessment of bisulfite conversion efficiency. After adapter ligation, bisulfite conversion was conducted using the EZ-96 DNA Methylation-Gold Kit (Zymo Research) according to the manufacturer’s protocol. PCR amplifications with 19 cycles were performed. Quality control of purified PCR products was performed on a 2200 TapeStation System (Agilent), and high-quality libraries were pooled and diversified with 15% PhiX. Single-end sequencing with 100-bp read length was conducted on a HiSeq 2500 sequencer (Illumina).
Quality assessment, data filtering, and mapping (RRBS)
In total, 106 individuals (48 wild-caught and 58 experimental fish) of balanced sex ratio were DNA-sequenced at an average of 19.8 ± 3.5 million reads for experimental fish and 11.4 ± 2.1 million reads for wild-caught fish (table S4). Demultiplexed Fastq files were quality-checked using FastQC version 0.11.5 (48) and MultiQC version 1.3 (44). Adapters were removed with cutadapt version 1.9.1 (49) using multiple adapter sequences (NNAGATCGGAAGAGCACAC, AGATCGGAAGAGCACAC, and ATCGGAAGAGCACAC) with a minimum overlap of 1 bp between adapter and read. This was necessary to remove primer dimers and avoid false methylation calls systematically caused by the RRBS end-repair step during library preparation, if the end-repair step adds artificial cytosines. Simultaneously, cutadapt was used to trim low-quality bases (-q 20) from the 3′ end and remove trimmed reads shorter than 10 bases. An air bubble during sequencing caused the bases 66 to 72 of 10 tiles of one lane (affecting 12 individuals) to have low-quality values, which were removed in a custom awk script. Two poor-quality individuals (a SYL and a NYN female) did not meet our strict quality requirements (e.g., ≥5 million reads; mapping efficiency, >52%) and showed biases in the proportion of bases per position compared to other individuals (plot in FastQC “per base sequence content”). Therefore, we excluded these two libraries from downstream analysis resulting in 15 instead of 16 individuals from SYL and NYN (Fig. 1). Bisulfite conversion efficiency was assessed from the spike-in controls (Cambridge Epigenetix) using the cegxQC software (50). Overall, conversion levels were 2.4 ± 1.8% conversion of methylated cytosines and 99.6 ± 0.5% conversion of unmethylated cytosines, which is in line with expected conversion rates (table S4). We used Bismark version 0.17.0 (51) to index the University of California Santa Cruz stickleback reference genome (Broad/gasAcu1) and to generate the bisulfite alignments with Bowtie2 version 2.3.3 at default settings. Bismark was also used to extract the methylation calls. Average mapping efficiency was 63.7 ± 2.4% (table S4).
Identification of differentially methylated sites
The methylation calls were analyzed in R version 3.4.1 (52) using the package methylKit version 1.3.8 (53). CpG loci were filtered for a minimum coverage of 10 reads per site. To account for potential PCR bias, we additionally excluded all sites in the 99.9th percentile of coverage. To improve the methylation estimates, we corrected for SNPs, which could have led to a wrong methylation call. The excluded positions were derived with custom-written Perl scripts from C-to-T and G-to-A SNPs with genotype quality of 20 and a minimum allele frequency of 0.005 (see above) from the 96 wild-caught individuals with a combination of custom-written Perl and R scripts using packages from methylKit (53) and GenomicRanges (54). After normalizing coverage values between samples, using normalizeCoverage implemented in methylKit, we excluded all sites that were present in fewer than nine individuals per treatment group from downstream analysis. As previously shown, sex-specific methylation affects <0.1% of CpG sites on autosomal chromosomes but >5% of CpGs on the sex chromosome (18). Therefore, to exclude a potential sex bias, we removed all CpG sites located on the sex chromosomes (chromosome 19), resulting in a high-quality dataset with 525,985 CpG sites. Last, by checking the first six principal components of the resulting principal components analysis and running an ANOVA on the filtered dataset, we confirmed the absence of an effect of sex on global methylation pattern (F124,1 = 2.611, P = 0.109). However, the principal components analysis revealed a bias in methylation pattern by families over all experimental groups. Therefore, to identify differentially methylated CpG sites (DMS) between treatment groups, we performed pairwise comparisons (table S5) fitting a logistic regression model per CpG site with calculateDiffMeth in methylKit using family as covariate for the experimental groups. A chi-square test was applied to assess significance levels of DMS, and P values were corrected to q values for multiple testing using the sliding linear model method (55). In addition, we accounted for multiple use of groups in pairwise comparisons and adjusted the α for the q value according to Bonferroni correction to 0.0125 (0.05/4). Ultimately, CpG sites were considered to be differentially methylated with a q < 0.0125 and a minimum weighted mean methylation difference of 15%. To ensure that the DMS obtained are not laboratory artifacts, we used calculateDiffMeth implemented in methylKit and compared the wild population from KIE to the experimental control group (KIE population from 20 PSU at 20 PSU). The resulting 11,828 DMS were excluded from the DMS obtained by the pairwise comparisons mentioned above (table S5). DMS were plotted across the genome for the comparison between KIE versus NYN (20 versus 6 PSU; blue fish) and KIE versus SYL (20 versus 33 PSU; yellow fish) using ggplot2 (56) and hypoimg (57) (fig. S1).
Assessment of inducibility and gene association of DMS
By comparing wild-caught individuals from the mid-salinity population (20 PSU; KIE) to the populations sampled at low (6 PSU; NYN) and high (33 PSU; SYL) salinity in the field, we obtained 1470 (KIE-NYN) and 1158 (KIE-SYL) pairwise pop-DMS. We first tested whether these pop-DMS distinguishing natural populations are inducible or stable at the respective salinity in the experiment. A pop-DMS was considered stable when the within- and the transgenerational acclimation groups did not significantly differ in methylation to the control group (q ≥ 0.0125). On the other hand, pop-DMS were considered inducible when at least one of the acclimation groups was differentially methylated compared to the control group (q < 0.0125; methylation difference, ≥15%). pop-DMS with a significant q value not exceeding the threshold of differential DNA methylation (15%) will be referred to as inconclusive hereafter. We used a randomization test to ensure that the number of inducible sites obtained did not occur by chance. To this end, we randomly sampled 1470 (KIE-NYN) and 1158 (KIE-SYL) pop-DMS from the complete dataset (1000 replicates). A chi-square test was used to assess whether our observed number of inducible, stable, and inconclusive sites differs from a random set of sites (averaged over replicates). Last, we tested whether the weighted mean methylation difference (meth.diff, in percentage) between wild populations matches the inducible methylation difference by subtracting the “meth.diff” in the experiment (exp) from the meth.diff between wild-caught populations (wild)
As we subtracted this difference from 100, values closer to 100 indicated higher similarity of experimentally inducible methylation with the postulated adaptive DNA methylation pattern in natural populations. By comparing the “δ.meth.diff” for within- and transgenerational acclimation using an ANOVA, we can assess whether there is a difference in inducibility of methylation to match patterns found in wild-caught populations. All analyses were run separately for decreased (6 PSU; KIE-NYN) and increased (33 PSU; KIE-SYL) salinity.
To detect potential functional associations of the observed changes in DNA methylation state, we classified the genomic region of a pop-DMS on the basis of their nearest TSS using annotateWithGeneParts and getAssociationWithTSS implemented in genomation version 1.4.2 (58). We distinguished between promoter (1500 bp upstream and 500 bp downstream of TSS), exon, intron, and intergenic regions. To be associated to a gene, the pop-DMS had to be either inside the gene or, if intergenic, not further than 10 kb away from the TSS. We excluded three pop-DMS that were on a different reference scaffold and then the gene they were associated to on the chrUn linkage group (that merges scaffolds into one large artificial chromosome). Using the genes with associated pop-DMS, we applied a conditional hypergeometric Gene Ontology (GO) term enrichment analysis (false discovery rate–corrected P ≤ 0.05) with the Ensembl stickleback annotation dataset “gaculeatus_gene_ensembl,” and all genes that were associated to any sequenced CpG site were used as universe. We identified overrepresented biological processes, molecular functions, and cellular components using the packages GOstats version 2.5 (59) and GSEABase version 1.46 (60) and corrected for multiple testing using the false discovery rate method implemented in goEnrichment version 1.0 (61) in R version 3.6 (52). Figures were produced using ggplot2 version 3.2 (56).
Estimation of DNA sequence–based genetic differentiation at differentially methylated sites
To evaluate the genetic differentiation up- and downstream (in sum, 10 kb) of the pop-DMS position, we calculated the mean FST values (≤60% missing data and depth, ≥5) from WGS data of the exact same individuals with vcftools version 0.1.15 (62). We hypothesized that inducible CpG positions matching the methylation difference expected from the profile of the wild populations are genetically more similar between the populations than sites that changed in the opposite direction. To test this one-tailed hypothesis, we applied a randomization test (with 10,000 bootstraps) on the mean FST difference between the two groups (expected and opposite)
We plotted the 10,000 delta mean FST values and calculated a P value by dividing the proportion of values smaller than the true difference by the number of bootstraps. Figures were produced using ggplot2 version 3.2 (56).
SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/12/eaaz1138/DC1
Fig. S1. Significant DMS throughout the genome for comparison between KIE versus NYN (20 versus 6 PSU; blue fish) and KIE versus SYL (20 versus 33 PSU; yellow fish).
Fig. S2. GO terms for biological processes, cellular components, and molecular functions under salinity increase (20 versus 33 PSU; yellow) and decrease (20 versus 6 PSU; blue) associated with pop-DMS.
Table S1. Relative distribution of DMS among genomic features.
Table S2A. Differentially methylated genes between populations from KIE (20 PSU) and NYN (6 PSU).
Table S2B. Differentially methylated genes between populations from KIE (20 PSU) and SYL (33 PSU).
Table S3A. Tukey post hoc test results for survival rate.
Table S3B. Tukey post hoc test results for SDL.
Table S3C. Tukey post hoc test results for HSI.
Table S3D. Tukey post hoc test results for total weight.
Table S4. Summary statistics for whole-genome resequencing of wild-caught sticklebacks.
Table S5A. Summary statistics for the RRBS of experimental fish.
Table S5B. Summary statistics for the RRBS of wild-caught fish.
Table S6. The number of DMS for each of the two pairwise population comparisons (pop-DMS).
This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
REFERENCES AND NOTES
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
S. A. Smallwood, G. Kelsey, Genome-wide analysis of DNA methylation in low cell numbers by reduced representation bisulfite sequencing, in Genomic Imprinting: Methods and Protocols, N. Engel, Ed. (Humana Press, 2012), pp. 187–197.
- ↵
- ↵
- ↵
CEGX Bioinfomatics Team, Cambridge Epigenetix (CEGX) (Babraham Research Campus, Cambridge, 2015).
- ↵
- ↵
R Core Team, R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2017).
- ↵
- ↵
- ↵
- ↵
H. Wickham, ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).
- ↵
- ↵
- ↵
- ↵
M. Morgan, S. Falcon, R. Gentleman, GSEABase: Gene set enrichment data structures and methods (2019).
- ↵
A. Hallab, goEnrichment: Helper functions to compute GO enrichment tests using GOstats and GSEABase (2015).
- ↵
- ↵
- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).