Contrasting DNA methylation responses of inbred fish lines to different rearing environments

ABSTRACT

Epigenetic mechanisms generate plastic phenotypes that can become locally adapted across environments. Disentangling genomic from epigenomic variation is challenging in sexual species due to genetic variation among individuals, but it is easier in self-fertilizing species. We analysed DNA methylation patterns of two highly inbred strains of a naturally self-fertilizing fish reared in two contrasting environments to investigate the obligatory (genotype-dependent), facilitated (partially depend on the genotype) or pure (genotype-independent) nature of the epigenetic variation. We found higher methylation differentiation between genotypes than between environments. Most methylation differences between environments common to both strains followed a pattern where the two genotypes (inbred lines) responded to the same environmental context with contrasting DNA methylation levels (facilitated epialleles). Our findings suggest that, at least in part, DNA methylation could depend on the dynamic interaction between the genotype and the environment, which could explain the plasticity of epigenetically mediated phenotypes.

KEYWORDS:

• Environmental enrichment
• epigenetic variation
• mangrove killifish
• RRBS
• genotype-by-environment interaction

Introduction

Epigenetic modifications are one of the potential molecular mechanisms to explain phenotypically plastic responses within genotypes [1,2]. This is because epigenetic markers can be altered by environmental variation and shape gene expression without changing nucleotide sequences [3], and ultimately affect phenotypic variation [1,4,5]. In phenotypic plasticity studies, the genome and epigenome are often experimentally confounded [6], and an implicit assumption is made that they react to environmental variation following similar norms of reaction. However, this may not necessarily be the case and genomes and epigenomes may respond differently to environmental change, thereby generating additional phenotypic variation [3,7].

To what extent epigenetic modifications act independently of genomic variation is key to understanding the potential role of epigenetics in evolution [3,7,8], as epigenetic variation completely under genetic control would not contribute any additional adaptive value [3]. Richards [9] classified epigenetic variation in obligatory, facilitated or pure epialleles, based on their degree of autonomy from the underlying genotype. Obligatory epialleles would be fully dependent on genetic variation and should show no variation across environmental change [10], whereas facilitated and pure epigenetic variation would differ in their degree of autonomy from the genotype (from partially dependent to independent) [10], acting as potential intermediaries between environmental conditions and genome responses.

Among the epigenetic modifications, DNA methylation is the best studied and plays an important role in the pre-transcriptional control of several biological processes, such as cell differentiation and genomic imprinting [11,12]. While correlations among DNA methylation patterns, environmental conditions and phenotypic traits have been widely investigated [8,13,14], the relative contributions of the genetic background and environmental variation to DNA methylation plasticity are still unclear [4,15]. Studies in humans and model organisms suggest that DNA methylation is influenced by the genotype, the environment and also by their interaction [15,16], but quantifying their relative influences is particularly challenging in natural populations with high levels of genetic variation [17]. In particular, there is little information on the basis of DNA methylation plasticity beyond model organisms [3], specially in teleost fishes, as most of the fish studies are done on zebrafish or salmonids [18].

The self-fertilizing hermaphroditic mangrove killifish (Kryptolebias marmoratus) [19] has naturally inbred lines [20] which inhabit mangroves that vary markedly depending on tidal conditions, ranging from temporary pools to mangrove leaf litter and crab burrows [21,22]. The species displays considerable plasticity in behaviour [23] and reproduction (mixed-mating with different degrees of self-fertilization and outcrossing) [24], both between and within self-fertilizing lines [25], and it has been suggested that regulation of gene expression through DNA methylation could play a role in its plastic response to environmental variation [26,27].

Inbred organisms provide a unique opportunity to detangle genetic from epigenetic variation [2,28]. Here, we investigated the relative roles of the genotype and the rearing environment (with or without physical enrichment) in DNA methylation plasticity of two genetically different and naturally inbred self-fertilizing lines of Krytolebias marmoratus. We hypothesized that if DNA methylation was mostly autonomous and shaped by environmental change, a higher number of different epialleles would be found between environments, regardless of the genetic background, than if DNA methylation was mostly under genetic control, where most of the epigenetic differences would occur between genotypes.

Results

We compared DNA methylation patterns in the brain of fish from two highly inbred mangrove killifish lines (DAN and R, originally collected in Belize mangroves but maintained under laboratory conditions for at least 20 generations) reared under physically enriched (with log and plants) or impoverished conditions (barren) for 10 months. By using Reduced Representation Bisulphite Sequencing (RRBS) we identified 5.5 million cytosine sites, of which 139.908 CpG sites fulfilled the minimum coverage requirement, representing 1.2% of the total number of cytosines of the mangrove killifish genome. This result is similar to recent RRBS studies in other fish (1% in rainbow trout [13]; 1.5–2% in guppies [29]).

The majority of cytosines surveyed mapped gene bodies (71.32%) or intergenic regions (19.10%), while only 2.54% were located on putative promoters. Linear models using the PCA scores for 1064 DMCs and 194 DMRs identified between genotypes and environments revealed that genotypes explained more of the variance for PC1 (54% of overall variation for both DMCs and DMRs) than environment (Figure S2; Table S1). When predictors were analysed individually, differences between genotypes also corresponded to a higher number of DMCs (817 vs 594, four DMCs shared) and DMRs (43 vs 17, no DMR shared) than differences between environments (Figure 1 and S2). Within genotypes, 357 and 3632 DMCs (25 and 373 DMRs) were identified between environments, for DAN and R, respectively. An additional analysis on three different subsets of six randomly selected R individuals (to match the number of DAN individuals) was carried out to assess possible biases due to differences in sample sizes. This analysis identified similar number of DMCs between lines, suggesting that the difference between lines was not due to sampling bias (Table S2). Unsupervised hierarchical clustering revealed distinctive methylation profiles between groups, except for comparisons between environments, where one and two individuals from the poor environment clustered with individuals from the enriched environment for DMCs and DMRs, respectively (Figure 1 and S3).

Figure 1. Heat map illustrating percentage of methylation for all differentially methylated cytosines (DMCs) identified (a) between genotypes, (b) between environments, (c) between environments for DAN strain, and (d) between environments for R strain (logistic regression q < 0.01 and |ΔM|>20%, and t.test p < 0.01) using unsupervised hierarchical clustering. Rows represent a unique CpG site and columns individual fish.

Twenty-five annotated DMCs and four DMRs between environments were shared across genotypes, potentially representing environmentally affected DMCs, independently of genetic background. Of these, based on the direction of methylation across environments, 22 out 25 DMCs were classified as potentially facilitated, with methylation scores following a genotype-specific pattern under similar environments (Table 1). This pattern was supported by the PCA results based on the DMCs methylation scores, which indicated different methylation profiles between environments (PC1 explaining 55.8% of variation), as well as genotypes (PC2 explaining 22.4% of variation) (Figure 2(a)). PC1 loadings were significantly influenced by the environment (t = 1.63, df = 1, p = 0.003) and the interaction between genotype and environment (t = −11.25, df = 1, p < 0.001), while PC2 loadings were only significantly influenced by the genotype (t = −1.64, df = 1, p < 0.001) (Table 2a). Methylation differences (with a lower threshold of 20%) for the facilitated DMCs ranged from 20.2% to 48.6% (Table 1). The potentially facilitated DMCs were mostly hypermethylated on enriched environments with respect to poor environments for DAN fish, while the opposite pattern was found for R fish (Table 1; Figure 2(c) and S4). The four DMRs between environments and shared by genotypes were also classified as facilitated, following the same methylation pattern found on the facilitated DMCs (Figure S5; Table S3)

Table 1. Methylation differences averaged (percentage) for differentially methylated cytosines (DMCs) between environments (poor, enriched), shared between genotypes (DAN, R) which overlap annotated genes (reference genome ASM164957v1, GCA_00164975.1). Epiallele classification (pure or facilitated) followed [9]. Positive and negative values represent increased and decreased methylation in enriched and poor environments, respectively. Q-value is the p-value adjusted by false discovery rate (FDR = 0.05).

Gene symbol	Entrez gene name	Epiallele classification	Meth diff DAN	Q-value	Meth diff R	Q-value
acvr2a	Activin A receptor type 2A	P	34.37	0.007	29.17	<0.001
col25a1	Collagen type XXV alpha 1 chain	F	43.61	0.005	−21.63	0.005
dmap1	DNA methyltransferase 1 associated protein 1	F	26.99	<0.001	−21.69	<0.001
foxp4	Forkhead box P4	F	22.50	<0.001	−37.56	<0.001
gpc5	Glypican 5	F	31.82	0.01	−32.06	<0.001
mipol1	Mirror-image polydactyly 1	F	35.85	<0.001	−30.59	<0.001
necab2	N-terminal EF-hand calcium binding protein 2	F	20.25	0.01	−23.56	<0.001
neo1	Neogenin 1	F	20.25	<0.001	−31.24	<0.001
nudcd1	NudC domain containing 1	F	39.76	<0.001	−21.05	<0.001
ramp3	Receptor activity-modifying protein 3-like	P	−27.12	0.037	−29.15	<0.001
ryr3	Ryanodine receptor 3	P	−30.48	0.003	−25.4	<0.001
sorcs2	Sortilin-related VPS10 domain containing receptor 2	F	36.81	0.008	−33.20	<0.001
trit1	tRNA isopentenyltransferase 1	F	20.38	<0.001	−29.70	<0.001
trmt44	tRNA methyltransferase 44	F	23.98	<0.001	−28.21	<0.001
ubald1	UBA like domain containing 1	F	36.76	0.019	−31.02	0.004
zeb2	Zinc finger E-box binding homeobox 2	F	31.83	<0.001	−28.35	<0.001
znf516	Zinc finger protein 516	F	31.51	<0.001	−31.60	<0.001
zranb3	Zinc finger RANBP2-type containing 3	F	41.41	<0.001	−33.79	<0.001
LOC108234847	Adhesion G protein-coupled receptor L3-like	F	48.57	<0.001	−39.35	<0.001
LOC108240988	Non-muscle caldesmon-like	F	32.49	0.005	−32.78	<0.001
LOC108243470	Protein-methionine sulfoxide oxidase mical2b-like	F	33.35	<0.001	−25.70	<0.001
LOC108243852	Receptor-type tyrosine-protein phosphatase N2-like	F	37.97	0.014	−30.28	<0.001
LOC108245430	Uncharacterized protein	F	42.94	<0.001	−41.49	<0.001
LOC108247402	Spectrin beta chain, non-erythrocytic 1-like	F	40.73	0.006	−27.58	<0.001
LOC108251479	Transcriptional regulator Myc-B-like	F	30.49	0.010	−25.13	<0.001

Table 2. Linear model of principal component scores for mangrove killifish epialleles shared between genotypes (R, DAN) and environments (poor, enriched).

	t-values	Prop. of variance (%)	df	p-value
(a) Facilitated epialleles
PC1 scores
Genotype	1.63	0.06	1	0.80
Environment	1.58	12.41	1	0.003
Genotype x Environment	−11.25	68.08	1	<0.001
PC2 scores
Genotype	−1.64	91.35	1	<0.001
Environment	8.28	0.21	1	0.49
Genotype x Environment	10.29	0.35	1	0.38
(b) Pure epialleles
PC1 scores
Genotype	2.28	13.18	1	0.008
Environment	−2.81	59.28	1	<0.001
Genotype x Environment	0.98	0.09	1	0.80
PC2 scores			1
Genotype	0.29	37.96	1	0.003
Environment	−1.82	1.14	1	0.56
Genotype x Environment	−2.97	0.67	1	0.65

Figure 2. Principal component analysis (PCA) and reaction norms of epialleles across genotypes and environments. PCAs were based on individual scores of methylation across either (a) facilitated or (b) pure annotated DMCs. Dark yellow for DAN individuals on enriched environments; light yellow for DAN genotype on poor environments; dark green for R individuals on enriched environments; light green for R genotype on poor environments. Each reaction norm represents the change on averaged methylation scores (in percentage) for (c) facilitated and (d) pure epialelles annotated DMCs across environments. Different colours represent the genotypes (yellow for DAN; green for R). Different shapes (d) represent different annotated DMCs. Epialelles were classified according to Richards (2006). Detailed information for each annotated DMCs methylation score across genotypes is available in Table 1.

Only three of the annotated DMCs within or neighbouring gene bodies were considered pure (Table 1; Figure S6). Average methylation differences for pure DMCs ranged from 25.4% to 34.37% (Table 1). The PCA only using pure DMCs showed a different pattern from the facilitated DMCs, with the PC1 separating environments explaining 72.38% of the variation, and the PC2 partially differentiating genotypes explaining 18.85% of the variation (Figure 2(b)). PC1 loadings were significantly affected by the environment (t = −2.81, df = 1, p < 0.001) and the genotype (t = −2.28, df = 1, p = 0.008), while PC2 loadings were only significantly influenced by the genotype (t = 0.29, df = 1, p = 0.003) (Table 2b, Figure 2(d) and S6).

Molecular network analysis revealed a highly connected network linked by genetic interactions and co-expression interactions, that was composed by 23 input-annotated DMCs (the uncharacterized LOC108245430 and ubald1 with no identified connections were removed) and 20 neighbouring genes (Figure S7). Centrality parameters, such as average degree (mean = 10.55; SD ± 5.89), closeness (mean = 0.53 ± 0.06), and radiality (mean = 0.77 ± 0.06) (Table S4), suggested that any alteration of the expression of the genes contained in the network might have major effects on genetic interactions and gene expression levels.

Twelve of the 15 most connected genes within the network (>10 connections) were input genes (i.e. genes affected by DMCs between environments and shared by genotypes). Gene ontology analysis showed that some of these genes are involved on important cellular and metabolic processes in zebrafish, such as regulation of transcription by RNA polymerase and gene expression (myc), RNA modification (trit1), intracellular calcium content (ryr3), and lipid metabolism (sorcs2), as well as pathways related to angiogenesis and stress response (ryr3 and myc) (Table S5).

Discussion

The potential adaptive role of epigenetically mediated plasticity depends on the relationship between the genome, the epigenome and the environment [6,7]. By using two naturally inbred strains of the mangrove killifish reared under contrasting environmental conditions, we have identified significant methylation differences among genotypes and environments, with different levels of autonomy from the genetic background.

Environmental enrichment affects the brain structure of fish [30–32], however few studies have investigated the molecular mechanisms underlying these changes [33,34] or whether the effects of enrichment vary across different genetic backgrounds. Kryptolebias marmoratus populations are composed by naturally inbred lines living in highly variable habitats [21,25], which display remarkable phenotypic variation (e.g. in reproductive output [24,35], behaviour [23,36,37] or sexual differentiation [35]), even under identical environmental conditions [25]. For example, the strains we used here were previously shown to display different sex-ratios in response to temperature variation [20], as well as differences in gene expression in response to parasitic infection [38], suggesting a potential combination of genetic and non-genetic mechanisms in mediating phenotypic variation [2]. Our results indicate that genotypes have an overriding influence on brain DNA methylation patterns and that their effect is greater than that caused by environmental enrichment. We only found a few DMCs that could be considered facilitated or pure epialleles, supporting the idea that environmentally induced autonomous DNA methylation may be limited [2]. Yet, the DNA methylation patterns of these putative independent epialleles indicated that DNA methylation outcomes could depend on specific combinations of the genotype and environmental conditions, although we cannot fully rule out the potential contribution of heritable epigenetic states independent of the genotypes [7] and/or brain cell heterogeneity.

The large differences in number of DMCs we found between lines could be explained by their genetic differences [39]. Studies in model organisms indicate that DNA methylation, and potentially other layers of chromatin organization, are strongly influenced by genomic variants [40–42]. For example, the spontaneous mutation in a gene related to methyltransferase1 activity increased by 40% the methylation differences among inbred lines of Arabidopsis thaliana [43], and in humans, 25% of variation in neonates’ methylomes can be explained by their genotype, while the remaining 75% is related to interactions between the genotype and maternal factors (i.e. smoking, age, intra-uterine environment) [15].

Most of the DMCs that we detected between environments and that were shared by both genotypes were located in gene bodies and were highly integrated within a gene network of genetic interactions and co-expression. Recent evidence in plants [44] indicates that gene body methylation can reduce erroneous transcription, and in oyster [45] and zebrafish [46] there seems to be a positive correlation between gene body methylation, gene expression and transcriptional regulation. Here, some of the genes affected by the DMCs found in gene bodies were related to the regulation of RNA polymerase activity and gene expression patterns (myc and trit1) [45,46], suggesting that these changes in methylation could be involved in biological and cellular processes.

DNA methylation is a good candidate for mediating phenotypic plasticity, given its responsiveness to environmental change, effects on downstream phenotypes, and transgenerational stability [3,4,47]. Our results suggest that, at least in part, DNA methylation patterns are influenced by a dynamic interaction between genotypes and the environment. Further research to investigate whether the patterns found here might influence transcription is warranted to assess the generality of our results, as this might provide a potential mechanistic explanation for the genotype-by-environment patterns often observed in phenotypically plastic responses [6,48].

Methods

We used hermaphrodite fish from two highly inbred strains (R and DAN) of Kryptolebias marmoratus originally collected from Belize [49,50] and kept under laboratory conditions (25–27°C, 16–18‰ ppm salinity under a 12 h light:12 h dark photoperiod), for at least 20 selfing generations [51]. The R (also called 50.91) strain was collected in Belize (Twin Cayes) in the early 1990s while the DAN (Dan06) strain was also collected from Belize in the early 2000s [49]. These selfing lines had previously shown different DNA methylation responses to environment (temperature) variation [20].

We compared brain methylation of fish reared under enriched and impoverished conditions, as previous studies had shown environmental enrichment can affect behavioural flexibility [52], brain size and cognition [53], and induce epigenetic modifications during early development [54]. We used two different habitats with different levels of environmental enrichment: 1) a physically enriched habitat, where individual fish were placed in contiguous rectangular tanks (9 cm depth x 12 cm width x 8 cm length) filled with 400 ml of brackish water and fitted with one artificial log (3 cm depth x 4 cm width x 4 cm length) and three artificial plants to simulate a complex habitat, and 2) a barren habitat (hereafter called poor) with the same tank conditions but without physical enrichment (Figure S1). Tanks were separated by opaque screens to prevent visual contact between individuals. For both strains, five initial lab-reared hermaphrodite parents of similar size (mean = 3.8 cm, SD = ±0.12) and age (mean = 417.3 days post-hatch, SD = ±13.4) were chosen. Eggs from these parents were maintained individually in circular plastic pots containing 100 ml of brackish water and checked daily (Figure S1). Upon hatching, individual alevins were randomly assigned to treatment tanks (enriched and poor), with one fish in each tank. Hatching success was 90%. The initial experimental set-up consisted of 29 R fish (18 in enriched habitat, 11 in poor habitat) and 21 DAN fish (10 in enriched, 11 in poor). Fish were maintained under standard laboratory conditions as above and fed three times a week with live brine shrimp (1 ml for the first two months post-hatch, and 2 ml thereafter). Fish were maintained in the experimental tanks for 10 months before being euthanized for brain methylation analysis. At 7 months post-hatch all fish in the experiment had laid at least one egg indicating that they were all sexually mature self-fertilizing hermaphrodites.

1. Genome-wide DNA methylation data

Fish were euthanized using tricaine methane-sulfonate (MS-222) following UK Home Office Schedule 1 and their brains kept in molecular biology grade ethanol (99%) for DNA extraction. Brain DNA was extracted from 22 individuals for epigenetic analysis (6 DANs: 3 from each environment; 16 Rs: 6 from poor, 10 from enriched environment) using Qiagen DNeasy Blood and tissue kit (Qiagen). Fish were genotyped for 23 microsatellites [55]. Genetic differences were identified between inbred lines (F_ST=1.00, Table S6), but not within lines. All individuals tested were homozygotes and identical within each line for all the markers analysed (Table S6).

Bisulphite converted genomic DNA libraries were prepared using Diagenode Premium Reduced Representation Bisulphite Sequencing (RRBS) kit according to manufacturer’s indications and sequenced on an Illumina NextSeq 500 platform using a 1x75pb single-end run, with PCR fully methylated and unmethylated spike controls added.

Quality assessment was performed using FastQC [56]. TrimGalore! [57] was used to trim low-quality base calls and adapters. Trimmed reads were aligned to the Kryptolebias marmoratus reference genome (ASM164957v1, GCA_00164975.1: source NCBI) prior in-silico bisulphite conversion using Bismark v0.17.0 [58], which was also used for cytosine methylation calls. Only methylation within CpG context [59] with a minimum coverage of 10 reads in each sample across the 22 individuals sequenced [60] was considered for subsequent analysis. Samples were divided into four experimental groups: ‘DAN enriched’, ‘DAN poor’, ‘R enriched’, ‘R poor’. Mapped reads were processed using SeqMonk [61]. After quality filtering, approximately 273 million reads were retained, averaging 12 million reads per sample. Of those, ~ 62.9% were uniquely mapped reads to the reference genome (Table S7). Overall bisulfite conversion was 99.6%.

1. Differentially methylated cytosines and regions

To identify differentially methylated cytosines (DMCs) across experimental groups, we used logistic regression on quantitated normalized data using α = 0.01 after multiple testing correction (Benjamini-Hochberg) and >20% minimal CpG methylation difference (|ΔM|), using R bridge in SeqMonk. We also performed t-tests across experimental group replicates, to generate a more conservative list of DMCs, only considering those shared by both statistical approaches. To identify differently methylated regions (DMRs), we performed a genome-wide unbiased DMR detection using tilling windows of 1000 bp on windows with at least five CpGs with ≥10 reads across all individuals.

We used the scores of methylation for DMCs and DMRs between genotypes and environments for principal component analysis (PCA) using ggfortify package [62] in R v. 3. 4. 3 (R Core Team 2014). To test for the effect of the genotype, environment and their interaction on the methylation scores, we used linear models with the scores for the first two PCA axis (>70% of the total variation) as a function of genotype, environment and their interaction. We then individually compared DMCs and DMRs between genotypes, followed by a comparison between environments. Subsequently, a comparison within each genotype between environments was carried out to identify potential environment-dependent DMCs and DMRs. From these comparisons, we identified annotated DMCs and DMRs shared between genotypes, which should represent commonly affected DMCs regardless of the genetic background.

We classified the DMCs and DMRs shared across genotypes between environments as facilitated, when displaying different directions of variation (non-parallel) on methylation scores across genotypes in the same environment (i.e. hypermethylated in an environment for one strain and hypomethylated in the other) or pure when displaying the same direction of variation (parallel) across genotypes and environments (i. e. hypermethylated or hypomethylated for both genotypes in the same environment) [9].

1. Molecular network analysis and centrality metrics

To identify potential functional implications of variation in DNA methylation for the annotated DMCs identified across genotypes between environments, we built a functional gene network using GeneMANIA [63]. To identify central genes [64] within the molecular network, we used NetworkAnalyzer [65] plugin into Cytoscape v. 3.7.1 [65]. Panther GO terms [66] was used to identify biological process and pathways for the most connected genes (>10 connections) within the network.

Authors’ contributions

SC, WMBF designed the experiment. WMBF, NB, DRB performed the experiment. WMBF, CGL, DRB analysed the data. WMBF, SC wrote the manuscript with participation of all authors.

Data accessibility

Sequences are accessible from https://www.ncbi.nlm.nih.gov/bioproject/506827.

Ethics

Work was carried out under Swansea University Animal Ethics Committee permit STU_BIOL_30484_110717192024_3.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

Work supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico, scholarship 233161/2014-7.