Effects of a parental exposure to diuron on Pacific oyster spat methylome.

Environmental epigenetic is an emerging field that studies the cause-effect relationship between environmental factors and heritable trait via an alteration in epigenetic marks. This field has received much attentions since the impact of environmental factors on different epigenetic marks have been shown to be associated with a broad range of phenotypic disorders in natural ecosystems. Chemical pollutants have been shown to affect immediate epigenetic information carriers of several aquatic species but the heritability of the chromatin marks and the consequences for long term adaptation remain open questions. In this work, we investigated the impact of the diuron herbicide on the DNA methylation pattern of spat from exposed Crassotrea gigas genitors. This oyster is one of the most important mollusk species produced worldwide and a key coastal economic resource in France. The whole genome bisulfite sequencing (WGBS, BS-Seq) was applied to obtain a methylome at single nucleotide resolution on DNA extracted from spat issued from diuron exposed genitors comparatively to control spat. We showed that the parental diuron exposure has an impact on the DNA methylation pattern of its progeny. Most of the differentially methylated regions occurred within coding sequences and we showed that this change in methylation level correlates with RNA level only in a very small group of genes. Although the DNA methylation profile is variable between individuals, we showed conserved DNA methylation patterns in response to parental diuron exposure. This relevant result opens perspectives for the setting of new markers based on epimutations as early indicators of marine pollutions.

Introduction

The role of environmental stressors in altering the physiological performance of aquatic wildlife is a fundamental topic that has recently attracted strong attention. Chronic environmental exposure to chemical compounds, including metals, pesticides, endocrine disrupting compounds can affect the behavior, and the physiological or developmental performance of aquatic organisms [1, 2]. Environmental epigenetics investigates the cause–effect relationships between specific environmental factors and the subsequent epigenetic modifications [3-7]. Chromatin structure is generally admitted as the epigenetic element of non-genetic inheritance system. Chromatin is a complex of DNA associated to proteins and epigenetic information carriers are the element which influences its structure. Non-coding RNA, methylation of cytosine residues and histone post translational modification are the main epigenetic information bearers known to affect chromatin structure. The epigenotype is mitotically and to some degree meiotically heritable, but unlike in the genotype, changes in the epigenotype are reversible. Therefore, the epigenetic information system (EIS, epigenetics) can provide heritable, novel phenotypes without variation in DNA sequence [8, 9]. The production of this heritable phenotypic variation may be essential for evolution [10] or may result in mal-adaptation of an organism experiencing changing environments as exemplified by several emerging disease disorders [11-13]. Various environmental contaminations have been shown to alter DNA methylation and histone patterns in animal species [14-17]. Currently, most of evidence of the association of altered epigenotypes and phenotypic disorders is provided from laboratory animal models and human epidemiological studies. Despite the extensive use of aquatic organisms in an ecotoxicological context, only a few studies provide information about the epigenetic mechanisms that may influence phenotypic variations in response to environmental factors in aquatic species [2, 18–22].

Amongst hazardous substances, diuron (N0-[3,4-dichlorophenyl]-N,N-dimethyl-urea) is a biologically active and quite persistent pollutant present in soil, water and sediments. This herbicide was used to eradicate weeds in non-agricultural areas (roadsides, railways, parks, etc.) [23]. Diuron is prohibited in France since 2008 as a phytosanitary substance and as a biocide since 2009 (Directive 2008/91/EC). But, since it is not well degraded naturally [24], it is still detected at high concentrations in French coastal and fresh waters [25]. Consequently, it is considered as a potential Persistent Organic Pollutants (POP), and is currently evaluated under the Stockholm Convention (http://chm.pops.int) [26]. Diuron’s ecotoxicological impacts have been studied on several aquatic species [27-29] including invertebrates of economic interest such as the Pacific oyster Crassostrea gigas. This species is reared in coastal areas (estuaries, lagoons and open sea) under the persistent threat of chemical pollution pressure. C. gigas is a successive and irregular protandrous hermaphroditic mollusc. Gonads are diffuse and non-permanent and their development is strongly influenced by seasonal seawater temperature variations upon which oysters release sperm and eggs into the seawater. Within 6 h, fertilized eggs develop into embryo and larvae in a planktonic stage lasting 2–3 weeks before entering metamorphosis leading to settlement to a solid substrate [30]. A few studies have investigated the gonadal and germline establishment in this species during development. However, following an orthologous to the vasa gene (called Oyvlg), Fabioux et al. [31] were able to follow primordial germ cells (PGCs) that would arise from 4 day blastomeres, hence suggesting that the PGCs specification mode in this species is preformation, i.e. maternally inherited [31-33]. Genomic and epigenomic resources for this species have recently become available [34, 35] and DNA methylation has been recently suggested to be involved in phenotypic variation in this organism [36]. The main goal of the present study was to determine how a parental exposure to the herbicide diuron can affect the DNA methylation pattern in C. gigas offspring. Some results were previously obtained demonstrating both direct and indirect effects of the parental herbicide exposure on oyster genome. A genotoxic effect has been detected through the evidence of primary DNA damages in both somatic, germinal and reproductive cells of diuron-exposed genitors [37, 38]. Moreover, the vertical transmission of damaged DNA was confirmed by the detection of chromosomal abnormalities in the offspring (DNA aneuploidy) [38]. By a transcriptional approach, we showed that exposure of oyster genitors to diuron at environmental concentrations affected gene expression including genes encoding enzymes involved in DNA methylation machinery [39]. In the offspring, a RNA-seq analysis demonstrated significant change in the transcriptional profile on F1 spat as well [40] that could explain developmental and growth impairment observed in F1 early life stages [37, 38].

In our study, we showed that a parental diuron exposure strongly affected spat DNA methylation pattern predominantly in intragenic regions and we discussed the possible consequences of this change on molecular phenotypes.

Results

Parental Diuron Exposure Causes Differential DNA Methylation in the Offspring

Whole Genome Bisulfite Sequencing was performed on three biological replicates for two investigated conditions (diuron-exposed (D) and seawater control (SWC). Sequencing of the samples yielded 620 166 410 100 bp Illumina paired-end reads, with an average of 103 361 068 paired-end reads by library (Supplementary Fig. S1). After quality control and alignment, an average of 36 965 778 (17.9%) reads by library was mapped (Supplementary Fig. S1).

On average, 16.6% of CpGs and 0.5% of CHH were found methylated, and the total methylation was 2.0% (5 mC of total DNA). There was no statistically significant difference in global methylation between D and SWC samples. The clustering analysis showed that the methylation profile obtained in exposure sample did not differ globally compared to the SWC as the different samples did not cluster according to their treatment (Fig. 1). For the three individual samples analyzed per group, parental diuron exposure did not globally modify the DNA methylation. Distance method used to produce the clustering is Pearson. We reasoned that if the epimutation rate would be higher in the diuron treated samples then the control samples should cluster together. Pearson correlation is also considered very sensitive to outliers. Strong differences in epimutation rates within diuron samples should produce such outliers, which is not the case. Consequently, differential methylation was restricted to specific regions. After visual inspection of methylation profiles we decided to use a window size of 500 bp for further analysis. Of the total of 1 124 112 500 bp windows of the genome in average 251 069 fulfilled the criterion of 10× coverage (Supplementary Fig. S1). However, in one sample, due to relatively poor coverage, only 57 819 windows were ≥10-fold covered. Within these 57 819 eligible windows, we found 1967 differentially methylated regions (DMRs) (3.4% of eligible windows): 1060 were hypermethylated and 907 hypomethylated in the D compared to SWC samples (Supplementary Table S2, sheet “MethDiff - Global results”).

Figure 1:

CpG methylation clustering. DNA methylation pattern clustering of the different samples is based on the Ward’s method and was performed with Bismark. Samples D3-17, D1-19, D3-12 correspond to spats from D samples TE1-14, TE3-19, TE1-11 correspond to spats from SWC

Target-Induced DMRs Occur Predominantly in Genes

Although the 1967 detected regions were statistically significantly differentially methylated between the D and the SWC individuals, it appeared that in many cases one individual among the three was less impacted. Therefore, we thoroughly examined the 1967 detected regions one by one to retain only those that are affected in the same way in the three biological replicates. This led to the selection of 236 DMRs (12% of DMRs, 0.4% of eligible windows). We therefore identified two types of epimutations: (i) environment-induced epimutations (in 0.4% of the epigenome, epimutation rate 4 × 10−3), and (ii) spontaneous epimutations (in 3% of the epigenome, epimutation rate 3 × 10−2). Following an earlier naming convention, we used the term “targeted epimutation” for type 1 epimutation, and “random epimutations” for type 2 [41]. Among the targeted DMRs, 121 were hypermethyled and 115 hypomethylated in D compared to SWC samples (Supplementary Table S2). The majority of these DMRs (73.2%) were found in genes, with 38.2% and 34.8% localized in exons and introns, respectively (Supplementary Table S2). Intergenic regions contain 30% DMRs among which 0.7% were found in putative promoter regions (only two DMRs). The distribution of hypermethylated and hypomethylated regions follows roughly the same pattern (Fig. 2, Supplementary Table S2). We found 170 DMRs located within genes of which 83 were hypermethylated and 87 hypomethylated (Supplementary Table S2, sheet “MethDiff—targeted epi hyper” and sheet “MethDiff—targeted epi hypo”). We further investigated the accurate location of the DMRs within the genes. 36.4% of the exon specific DMRs occurred in the last exon genes. 25.5% of the intron specific DMRs occurred in the last intron genes (Fig. 2). This overrepresentation of the DMRs was not observed in the first exon of genes if we counted the number of DMRs from the first exon position (data not shown). Therefore, we observed that DMRs preferentially occurred at the end of the coding sequences.

Figure 2:

Genomic feature of the target DMRs. The genomic features of the target DMRs were determined using as a guide the gene annotation obtained from previous transcriptomic data [40]. Hypo/hypermethylation refers to methylation level in D as compared to SWC. A and B represent genomic features for hypermethylated and hypomethylated regions, respectively. C and D represent the number of DMRs per exon or intron, respectively, counting from the last exon or intron position, dark and pale grey represents the hypermethylation and hypomethylation, respectively

Interestingly, we found several genes strongly affected by the parental diuron exposure (Fig. 3 and Supplementary Fig. S2). For these candidates, several 500 bp DMRs accumulates within the coding sequence. XLOC_005912 (CGI_10018082) encodes an uncharacterized protein with two conserved domains, one found in chromosome partitioning proteins and another one found in mitochondria eating proteins. It contained 63 intragenic CpG sites in a 1.437 kb window which were strongly methylated in the three SWC biological replicates. All these sites were poorly methylated in the three samples from the D individuals. XLOC_030549 (CGI_10011718) encodes a serine protease with domains found in fibrinolytic proteins. 16 strongly methylated CpG sites were found in a 1.132 kb window upstream of the transcriptional start site in the SWC group. These sites were unmethylated in the D group. The scaffold309 from position 131250 to position 137240 contains 67 CpG sites highly methylated in the individuals from D group compared to SWC group. This region contained four transposon-like elements (XLOC_019653, XLOC_019654, XLOC_019655, XLOC _019656) for which transcripts had been detected in previous RNA-Seq experiment [40]. In summary, among the 1967 differentially methylated regions, we found 12% which were targeted epimutations. These epimutations mostly occured within coding sequences and preferentially at the end of genes.

Figure 3:

Example of strongly impacted regions. Extract of the C. gigas jbrowse interface which allows visualization of transcriptomic and DNA methylation data for three regions: (A) Scaffold 1255 from 364500 to 370000 containing XLOC_005912 (CGI_10018082), (B) Scaffold 43664 from 222500 to 235000 containing XLOC_030549 (CGI_10011718), and (C) Scaffold 309 from 130000 to 137500 containing four unique exon genes (XLOC_019653, XLOC_019654, XLOC_019655, XLOC_019656). Upper panel represents the annotation of genes in the current genome assembly V9 (CGI number), and the CDS according to previous transcriptomic data [40]. Bottom panel represents the methylated cytosine position. X-axis represents the position of each methylated cytosine. Y axis (bars) represents the % of methylation observed for each methylated cytosine (scale is 0–100). The six tracks represent the three biological replicates of individual spats from genitors treated with diuron or control (from top to bottom)

About 80% of Environmentally Induced DMRs Are Specifically Targeted

As outlined above we identified 236 targeted DMRs using an epigenome-wide screen in single nucleotide resolution in three individuals. However, roughly 88% or all detected DMRs were random epimutations. It was therefore conceivable that all DMRs were actually random and just by chance, overlapping existing DMRs in our three replicates that were mistakenly assigned to be environmentally induced. We therefore randomly sampled 10 candidate DMRs among those that showed the most pronounced differences in DNA methylation and performed DNA methylation analysis on offspring spat of 10 individuals of each group (D, SC, SWC) by targeted bisulfite sequencing (T-BS-Seq) (Table 1, Supplementary Fig. S2). For 8 out of 10 samples, we observed a correlation between the BS-Seq and target analysis (Supplementary Fig. S2). In one loci (Scaffold42366: 107142.107725, Supplementary Fig. S2) we found good correlation between BS-Seq and T-BS-Seq for the three individuals studied by BS-Seq but not for the others. In one candidate (Scaffold1720:210182.210722), we did not confirm the BS-Seq results. Moreover, we observed several scenarios of diuron/solvent effect (Table 1 and Supplementary Fig. S2): a diuron specific effect was observed in a unique region. Otherwise, we observed a solvent specific effect in another region, a synergistic/same/antagonist effect was observed in two, three and one regions, respectively. And finally, two regions displayed really complex profile where the diuron/solvent mixture effect relied on each CpG site individually. Finally, we barely found the same DNA methylation profile between individuals whatever the environmental conditions (Supplementary Fig. S2). Despite this individual variation, we found two consistent loci for which DNA methylation pattern was really conserved between individuals. The diuron/solvent suppressed the methylation of all the CpG sites in all the tested individuals in the loci XLOC_005912/CGI_10018082 present on scaffold1255 (365000–370000), whereas this mixture induced an overall increase in DNA methylation on several CpG sites of a region spanning 3 kb of the scaffold33832 (12000–15000). These consistent results within the 10 tested individuals led us to calculate a positive predictive value (PPV). The PPV tests the performance of a diagnostic based on the proportion of positive and negative results. It aims at checking if the two candidate regions can be used as potential environmental contamination marker based on DNA methylation profile. A PPV of 85.7% was obtained when taking into account the combined results of these two loci (Supplementary Fig. S3).

Table 1:

DNA methylation analysis on target loci

Targeted loci	Diuron/solvent effect
Scaffold1154:309231.309691	Diuron specific effect
Scaffold433: 896527.897085	Solvent specific effect
Scaffold33832: 13205.13768	Solvent effect, further increased with diuron
Scaffold1255: 75.750
Scaffold41174: 81270.81829	Same effect for solvent and diuron
Scaffold42366: 108172.108673
Scaffold43170: 116309.116863
Scaffol37178: 4921.5527	Antagonist effect between solvent and diuron
Scaffold1720: 210182.210722	Diuron or solvent effect rely on CpG sites
Scaffold42366: 107142.107725

Ten loci were selected as indicated in column 1 and targeted by PCR after DNA bisulfite conversion. Ten individuals were analyzed per condition (SWC, SC, D). Detailed results are given in Supplementary Fig. S2.

No Indication of Genetic Selection Following Parental Diuron Exposure

Diuron exposure could have led to selection of genetically encoded phenotypes and the observed DNA methylation (epiallele frequency) differences could be an expression of such underlying genetic (allele frequency [AF]) difference. We used the genetic variants obtained with MethylExtract to search for signs of selection on a genome wide level. We used a total of 11 425 759 SNPs (control) and 9 517 133 SNPs (diuron) to calculate Tajima’s D in 200 kb sliding windows. Genome wide scan’s of Tajima’s D in sliding windows can identify regions whose D greatly deviates from the bulk of the empirical distribution and such windows are reported as significant. Frequency distribution of Tajima’s D indicates that most 200 kb windows have a value of around 0 and that there are no major differences in D distribution between SWC and D conditions. For a high-resolution analysis, we also cloned and sequenced a 2.198 kb amplified cDNA fragment from the locus XLOC_005912/CGI_10018082 in which DNA methylation is strongly affected by the diuron parental exposure. 248 SNPs and 14 Indels were detected in this 2.198 kb fragment of XLOC_005912/CGI_10018082 cDNA. The clustering analysis performed groups according to individuals rather than condition (Supplementary Fig. S4). While we could not formally exclude that minor changes in allele frequency occurred during parental diuron exposure, it is sure that there was no strong and systematic AF distortion.

Relationship between DNA Methylation and RNA Abundance

After having firmly established that parental diuron exposure led to either random or targeted changes in DNA methylation profiles we wondered if such changes could affect the level of RNA abundance. To establish the regions of methylation enrichment around genes we produced average methylation profiles for all 42 463 genes (Fig. 4). As expected, gene bodies were highly enriched for methylation and no difference could be detected between SWC and D offspring. Visual inspection of genome-wide methylation profiles confirmed the mosaic-style methylation with large blocks of high methylation separated by blocks with no or weak methylation. We established the following criteria for both methylation types: regions with at least 5 CpGs that are all at least 70% methylated, and having a maximum distance of 2 kb between them (read coverage ≥ 10) were considered as “high methylation region” (HMR). All other regions were considered as “low methylation region” (LMR) (Supplementary Fig. S1). HMRs contain 5 655 genes (13%) whereas LMRs contain 36 808 genes (87%). We then used our earlier RNA-Seq data [40] to evaluate the relationship between gene body DNA methylation and RNA abundance. For genes located in HMRs, we observed a weak correlation (Pearson coefficient 0.4) between degree of methylation and RNA abundance (Fig. 5). No correlation was found for genes located in the LMRs. For genes in the high methylation compartment, it was therefore in principle conceivable that changes in methylation could translate into changes in RNA level (or vice versa). Among the top 25 483 differentially expressed genes (abs(log2) fold change ≥2) only 80 were located in HMRs. Within these genes we observed two populations: (i) those that showed only moderate changes in RNA abundance and this despite strong changes in DNA methylation, and (ii) a smaller fraction where increase in DNA methylation correlated with increase in RNA amounts (Fig. 6). Amongst the above mentioned conserved 170 DMGs, high methylation compartment contains 84 (49%) (54 hypomethylated after diuron, and 30 hypermethylated). However, their abs(log2) expression differences were <2. In summary, it appeared that DNA methylation changes correlated directly with RNA level only in a very small group of genes, and only for genes located in the high methylation compartment of the (epi)genome. For this analysis, we used pooled BS-Seq data and compared them with pooled RNA-Seq data. A potential caveat of this approach is that potential individual differences in RNA abundance and DNA methylation were leveled out by absence of such differences in the replicates. Therefore, we performed qPCR experiments targeting arbitrarily selected candidate genes in DMRs using cDNA from the same individuals used for the BS-Seq experiment. We chose different contexts of DMR position: promoter, within genes, 3′ distant regulatory region (Table 2). Among the 12 tested genes, only two (XLOC_005912/CGI_10018082 and XLOC_032365/CGI_10010585) showed a significant difference in their expression levels between control and parental diuron exposure. XLOC_005912 was 1.80 times more expressed in the D than in the SWC sample. This gene was hypomethylated in the D sample all along its coding sequence (Table 2 and Fig. 3). Reversely, XLOC_032365 was 1.66 times less expressed in the D than in the SWC sample although it was hypomethylated in the D sample in an intronic region (Table 2). The other candidates did not show significant differences in their expression profiles although DNA methylation pattern was sometimes strongly affected by the parental diuron exposure (Fig. 3).

Figure 4:

Average DNA methylation for transcripts. The average DNA methylation pattern around genes in offspring from diuron exposed parents (A) and from sea water control parents (B) was generated using Deeptool based on 67 754 transcripts produced earlier by RNA-Seq [40]. X axis represents the gene context; gene body was scaled to 2 kb. Y axis represents the average percentage of DNA methylation

Figure 5:

Correlation between DNA methylation and gene expression. Correlation was performed using the HTSeq counts from previous RNA-Seq experiment [40] and average DNA methylation level per gene from the three biological replicates of each condition

Figure 6:

Correlation between the change in methylation level and change in expression level. The fold change in expression observed between the diuron and sea water control conditions (output from DESeq analysis [40]) was compared with the fold change in DNA methylation between the D and SWC conditions

Table 2:

Effect of DNA methylation on gene transcription

CGI number	XLOC number	Annotation	DMR positions	Context	Fold change
CGI_10018082	XLOC_005912	Centromeric protein	Scaffold 1255	Intronic and exonic	1.80
			366001–367000	Hypomethylation
CGI_10026162	XLOC_004001	Thap domain-containing protein 9-like	Scaffold 1154	Exonic and intronic	No significant difference
			309001–311000	Hypermethylation
	XLOC_020586	Myb-related transcription factor, partner of profilin	Scaffold 33832	Exonic and 3′ distant reg. region	No significant difference
			12001–12500	Hypermethylation
CGI_10009103	XLOC_029410	Hypothetical protein	Scaffold 43170	Exonic	No significant difference
			116501–117000	Hypomethylation
	XLOC_030549	Hepatocyte growth factor-like protein	Scaffold 43664	Promoter	No significant difference
			232001–232500	Hypermethylation
CGI_10011372	XLOC_017268	Reverse transcriptase	Scaffold 206	Exonic	No significant difference
			81001–81500	Hypomethylation
CGI_10010582	XLOC_032353	Synaptic vesicle glycoprotein 2c-like	Scaffold 473	Exonic	No significant difference
			121501–122000	Hypomethylation
CGI_10010585	XLOC_032365	28s ribosomal protein mitochondrial-like	Scaffold 473	Intronic	0.60
			130001–130500	Hypomethylation
CGI_10017761	XLOC_037874	Beta-arrestin1	Scaffold 713	Exonic (3′UTR)	No significant difference
			415001–415500	Hypomethylation
	XLOC_019656	Retrovirus polyprotein	Scaffold 309	Exonic	No significant difference
			136001–136500	Hypomethylation
CGI_10026723	XLOC_019464	Paired box pox meso protein	Scaffold 301	Exonic	No significant difference
			516001–516500	Hypermethylation
CGI_10023654	XLOC_034086	Solute carrier family 13 member 2	Scaffold 54	Intronic	No significant difference
			79501–80000	Hypermethylated

Genes have been chosen according to different context of DMR position. DMRs context are indicated in the “context” column. The accurate position on C. gigas genome assembly V9 is indicated in column “DMRs position.” Hypo/hypermethylation and difference in expression refers to D compared to SWC. Fold changes are indicated only in the case of significant P value (P < 0.01).

Since it has been suggested that intragenic DNA methylation plays a role in splice variant regulation [42], we investigated the occurrence of different splice variants between the D and SWC group. We amplified by conventional PCR long fragments that covered almost the full CDS for three candidate genes (XLOC_005912/CGI_10018082, XLOC_032353/CGI_10010582 and XLOC_032353/CGI_10010585) in 15 individuals per condition including the individuals which had been used for the BS-Seq experiment. The size of amplified products was checked on agarose gels and we observed a unique fragment of the expected size in both conditions for each targeted locus (data not shown), suggesting that differential methylation did not systematically lead to splice variants. Taken together, our data indicated that although the diuron exposure strongly affected the DNA methylation profile of spat from exposed genitors, differential methylation had generally no direct impact on RNA abundance and splicing.

Discussion

A growing body of evidence shows that environmental stressors can be responsible for changes in epigenetic marks, which may result in heritable phenotypic outcomes and changes in life history trait of living species [43-45]. The increasing data regarding the influence of the environment on chromatin structure has led to the emergence of the “environmental epigenetics” discipline [7, 19, 46, 47]. Recently, this field has received many attentions since DNA methylation and histone modifications have been shown to be associated with a broad range of disease outcomes, including cancer, metabolic disorders, allergies, asthma, and reproductive disabilities [16, 48–50]. Although the majority of the efforts are oriented to the context of human health, growing evidences suggest that epigenetic based events play also a major role in natural ecosystems [18, 51]. Substantial efforts have been devoted to the study of multigenerational impacts of the climate change (ocean acidification and temperature increase) [52, 53] and chemical pollution on several marine species [2], including mollusks [54, 55]. Diverse examples illustrate how parental exposure could benefit to the offspring which exhibit traits with increased fitness in the environmental conditions experienced by their parents [54, 55]. Other studies have shown how chemical pollutants such as endocrine disrupting compounds can have a negative impact on different life history traits such as fertilization success, egg viability, reproductive behaviors, hatch success and development [2, 56]. Although, there are evidences for an impact of the environment, including chemical pollutants, on DNA methylation content of several aquatic species [57-59], no studies have investigated the multigenerational or intergenerational effect of pollutants on epigenetic marks of the offspring of aquatic relevant species. In this work, we investigated for the first time the impact of a parental exposure to the herbicide diuron on epigenetic marks of C. gigas. The effect of this herbicide on different C. gigas life history traits has been demonstrated several times on exposed individuals but also on the progeny of exposed parents [37, 38, 60–62]. In the present study, we showed by BS-seq analysis on a restricted number of samples the existence of 1967 differentially methylated regions (DMRs) after a parental diuron exposure. These DMRs could be produced by the parental stress, but some of these differences could be the product of an interindividual variability of DNA methylation in the spats used. Among the 1967 DMRs, 88% are spontaneous epimutations and therefore the results of interindividual variability whereas 12% of these methylation differences are target-induced epimutations and therefore those regions are the results of the diuron-solvent mixture effect on the DNA methylation pattern. We then established the CpG methylation profile of 10 regions in a more subsequent number of samples by conventional PCR target analysis after DNA conversion with a bisulfite treatment. This analysis reveals the complexity of the diuron/solvent mixture effect on C. gigas methylome as we barely observed the same consequences depending on the considered loci, CpG within each loci and individuals. A diuron specific effect was observed in an unique case and we rather showed that the differential CpG methylation pattern were the result of different scenarios, including solvent effect only and a synergistic or antagonist effect of the diuron/solvent combination. Although we observed a strong individual variation in response to the exposure, two regions were consistently affected in the 10 tested individuals by the solvent/diuron combination. Using the combination of DNA methylation in both these markers, we calculated a positive predictive value which suggests that these two regions could indicate pollution in marine environment by a confidence level of 85.7% and could be used as potential biomarkers for marine risk assessment. The use of DNA methylation patterns as an early indicator of general disease disorder is already applied successfully in cancer diagnostics [63-65]. Further application as biomarkers in the evaluation of the carcinogenic potential of environmental factors has also been suggested [66]. The same kind of application in an ecologically relevant context to detect contamination in the environment on sentinel species has been already discussed but never applied [19, 51, 67]. In this context, these two regions represent relevant candidates which could be tested for the biomonitoring of pollution in the marine environment although further appropriate validations are required. Especially, it will be important to study whether the impact of parental diuron exposure on the DNA methylation profile of its progeny mirrors a general stress effect or corresponds to specific effect of this class of chemical pollutants.

Most of the studies conducted so far in the field of environmental epigenetic have been centered on DNA methylation. This covalent modification is conserved in most major eukaryotic groups [68-70], although it has been lost in certain model species [69, 71, 72]. The levels, patterns and function of DNA methylation vary among species. In invertebrates, it occurs in many taxa [69, 73–79] whose genomes are characterized by interspaced regions of methylated and unmethylated DNA. Intragenic (exon and intron) methylation is a general feature of invertebrate animals while methylation of repetitive elements and intergenic regions occurs only at moderate level [80]. The role of intragenic DNA methylation in invertebrate species remains an open question. It is more abundant in the 5′ region of open reading frame and associated with histone marks indicative of open chromatin structure in insect species [74, 81]. Our work confirmed previous data on the DNA methylation pattern in C. gigas [35, 77]. We showed that 2% of total cytosine was methylated predominantly in the CpG context, DNA methylation was organized as a mosaic feature and was higher in the coding region. Furthermore, our work highlighted that this pattern was sensitive to environmental condition and could be involved in intergenerational effects. The majority of DMRs that we identified was essentially found in intragenic regions. We investigated the consequences of these DMRs which occurred within genes on their corresponding transcripts. Some studies have shown that intragenic methylation is likely to occur in genes associated with stable transcription [74, 80] while other studies have suggested that DNA methylation may regulate gene expression by altering the landscape of splice variants [73, 82]. In our case, we observed that DNA methylation changes correlated with RNA level only in a very small group of genes, located in the high methylated regions. The target expression analysis confirmed the absence of obvious impact on gene expression and splice variant generation, although, some of these candidates displayed a striking difference in their CpG methylation profile and in one case a strong difference was observed in a promoter region. This lack of obvious phenotype raised the question of the function of DNA methylation for molecular phenotypes, if not playing a role in controlling gene expression as clearly demonstrated in flowering plants and mammals [68]. Another striking result of our study was the severe impact of the diuron parental exposure on several unique exon genes that encode transposon relic. These mobile genetic elements seemed to be active since we could detect their transcripts. However, there were no differences in their level of expression between the SWC and the D groups which further questions the role of CpG methylation for transposon silencing in C. gigas. Previous work has shown that DNA methylation only targeted some classes of repetitive elements in C. gigas [77] and therefore mobile elements were not a preferential target of DNA methylation as observed in plants and mammals for which their role as a genome defense has been largely demonstrated [68]. Consequently, our study confirmed previous observations that DNA methylation in invertebrates do not follow the general paradigm concerning its role for gene silencing [83-85]. We suggest that the admitted role of DNA methylation for heterochromatization of DNA must be revisited and its function in controlling gene expression in invertebrates, if there is any, may be restricted to really conserved pathway such as those encountered for developmental and key signaling pathways [86-89]. Our analysis further highlighted that a bias exists in the position of the DMRs which preferentially occurred in the 3′ end region of the genes. Given that little information exists on the possible role of invertebrate intragenic DNA methylation, the meaning of such a bias certainly deserves further attention. An important issue to address is whether the observed effect is linked to the genotoxic effects of diuron observed in parallel during this experiment or whether it is a general feature that is observed whatever the stress experienced. As a matter of fact, some studies pointed out a possible interaction between genotoxicity and DNA methylation mechanisms [90].

In conclusion, we have investigated the impact of a parental diuron exposure on the DNA methylation profile of C. gigas. We showed that a parental chemical stress effect can have an impact on the DNA methylation pattern of the progeny. This is the first evidence for the implication of epigenetic marks for a parental effect in aquatic species. Our work showed that epigenetic changes might represent an important pathway by which environmental factors influence phenotypic outcomes, within individuals and across generations in invertebrate aquatic species of economic and ecological importance.

Methodology

Genitor Origin, Parental Diuron Exposure, and Production of the Next Generation

The biological material used for the present study was obtained at the experimental hatchery of Ifremer La Tremblade in the frame of the GIMEPEC project funded by the French National Agency for Research (ANR-CESA-01601) [37]. Briefly, wild adult oysters from Marennes-Oléron Bay (France) were transferred from the field to the hatchery. Seawater was pumped directly from the Seudre river Estuary, filtered through a sand filter (40 μm) and passed through UV rays before draining the tanks in a continuous oxygenated flow system. Oysters were fed daily with a mixture of four marine microalgae and the water temperature was maintained at 8°C (±1°C) throughout the 1 month acclimatization period. In order to initiate the oyster gonad maturation process, the temperature was then raised by two degrees per day for 1 week, to reach 19.8°C (±0.3°C) at the start of the experiment. A total of 720 individuals were divided into three experimental groups (three replicate tanks per group): a diuron-exposed group (D), a solvent-exposed group as diuron was prepared in 0.005% acetonitrile (SC) and a seawater control group (SWC). At the start and the mid-course of gametogenesis, the oysters were exposed during two 7-day periods to pulses of 0.2 and 0.3 µg/l of diuron, respectively, and of 0.005% (50 µl/l) acetonitrile for the SC group. The effective diuron concentrations reached into the assay tanks were measured thanks to the use of Polar Organic Chemical Integrative Samplers (POCIS) [37]. This exposure scenario is expected to be environmentally realistic with an exposure time which mimics rain events at diuron concentrations close to those detected in coastal waters [25, 91]. In France, Diuron has been detected at concentrations up to 0.25 µg l−1 [92]. At time of spawning, male and female gametes from each genitor lots (“SWC”, “SC” and “D”) were recovered and used to produce the next generation as previously described [37].

Genomic Isolation

Individual whole tissue samples of spat oysters previously frozen in liquid nitrogen were homogenized in RL1 lysis buffer supplied by the NucleoSpin 96 RNA Tissue Core kit from Macherey-Nagel. Genomic DNA (gDNA) was isolated individually as previously described [93]. Briefly, 100 µl of oyster tissue were incubated overnight at 55 °C in 500 µl of DNA extraction buffer (100 mM NaCl, 10 mM Tris–HCl pH 8, 25 mM EDTA pH 8, 0.5% SDS, 0.1 mg/ml protease K). Following a step of phenol/chloroform extraction, gDNA samples were precipitated by addition of cold ethanol and treated with 50 µg/ml RNase A (Invitrogen) for 30 min at 37 °C to eliminate contaminating RNA.

Genome Wide Methylome Analysis

gDNA (3.02–5.44 µg) was sent to GATC (Germany) for bisulfite conversion, BS-seq paired-end library construction and sequencing. Bisulfite conversion rate was 99.5%. Three libraries per oyster group (SWC and D) were constructed as biological triplicates. One library corresponds to whole DNA sample from a single individual. Individual D3-17, D1-19, and D3-12 were chosen for the D group and individuals TE1-14, TE3-19, and TE1-11 were chosen for the SWC group. These individuals are from the same study for which we previously generated a RNA-Seq experiment except that in this RNA-Seq study, reverse transcription has been realized on pool of RNA from 10 individuals including the one analyzed here [40]. The workflow applied for BS-Seq analysis is shown in Supplementary Fig. S1. Fastq reads of six libraries (three per condition) were filtered separately based on quality phred-score. Reads with <95% of nucleotides with a quality phred-score equal or higher than 26 were discarded with the Filter by quality tool. Interlace of forward and reverse libraries per sample were performed to discard unpaired reads. De-interlace step was performed per sample to separate forward and reverse libraries previous to mapping. The reference genome, CgigasV9 [34], was modified for the subsequent mapping step by discarding contigs ≤ 5 Kb. Bismark [94] was used to map the filtered reads to the modified genome with zero mismatch allowed in the 100 bp reads. Bismark sorted alignment output (BAM file) were converted to SAM file on a local Galaxy instance. Coverage was analyzed using “bedtools makewindows” [95] and “samtools bedcov” [96]. The methylkit software (R packages) [97] was used for comparative statistical analysis between the samples using the SAM files as input. Criteria used to determine differentially methylated regions (DMRs) of 500 bp were as follow: a minimum coverage of 10× in all replicates, a minimum of 25% of methylation difference between conditions with a corresponding adjusted P value (FDR) of 0.01. Methylkit defines methylation difference as absolute value of methylation percentage change between test and control. Since there is no data in oyster about the degree of methylation change that is necessary to obtain a biological impact we tentatively used the 25% default value. MethylExtract 1.8.1 [98] was used to generate Methylation Wiggle and BED files that were displayed in IGV [99] for visual inspection of the data genomic regions (promoter, exon, intron or intergenic region) affected by parental exposure to diuron. MethylExtract served also to generate genetic variant calls in VCF format. VCF files were analyzed with VCFTools v0.1.12b [100]. Average methylation profiles over transcribed regions were generated with DeepTools [101] on the public Galaxy instance (https://usegalaxy.org/). We used our earlier produced RNA-Seq based annotation of the genome [40]. Correlation analysis between DNA methylation and gene expression was performed using the Pearson and Apos correlation tool under the public galaxy instance (https://usegalaxy.org/). Average HTSeq count per gene [40] was used for the gene expression level. Average DNA methylation profile per gene for the three biological replicates was generated. The fold change in expression observed between the diuron and seawater control conditions was the result of the DESeq analysis [40]. The fold change in DNA methylation between the diuron and seawater control conditions was the result of the Bismark analysis. For the bioinformatic pipeline, scripts are available upon request.

Visual inspection of the methylation profile of each replicates in IGV software was performed in order to identify DMRs for which all biological replicates were affected in the same way. We determined the genomic features of these target DMRs as being localized either within coding region (exon or intron), promoter or intergenic region using the gene annotation obtained with previous RNA-seq data [40]. A 500 bp DMRs which was encompassing two genomic features were counted as being located in the two contexts. Promoter regions were considered to be the 1 kb region upstream of the first exon of an annotated gene. For accurate exonic position within the genes, we defined 13 possible positions: position in unique exon, in last exon, in last exon −1 to last exon −10 and others which concerns genes with more than 10 exons. The same definition was used to localize the accurate intronic position within the genes.

DNA Methylation of Target Gene Analysis

Bisulfite treatment was done as described before [102] (http://ihpe.univ-perp.fr/methods/methods/bisulphite_engl_mc.html) on gDNA extracted from SWC, SC and D groups. Ten individuals were analyzed for each experimental condition. Briefly, 300 ng of DNA was denaturated with 3 M NaOH, treated with a solution of sodium-bisulfite and hydroquinone at pH 5 in the dark for 4 h at 55 °C; desalted (Amicon Ultra column, UFC510024 Millipore), desulfonated by adding 350 µl of 0.1 M NaOH to the DNA in the Microcon column, and finally dissolved in 10 mM Tris/Cl pH 8. The genomic sequence of 10 DMRs was targeted by nested PCR on 10 individuals per group (SWC, SC, D). Primer pairs (Supplementary Table S1) were designed using MethPrimer [103]. PCR amplification was performed in 50 µl using 1.25 U of the M0495L-Hot start Taq DNA polymerase (Biolabs), 1× Hot start Taq buffer, 0.2 µM of primers, 200 µM of dNTP. PCR cycles were:94 °C for 2 min, five cycles of 94 °C for 1 min, 50 °C for 2 min, and 72 °C for 3 min; followed by 25 cycles of 94 °C for 30 s, 50 °C for 2 min, and 72 °C for 1 min 30 s and finally 72 °C for 10 min on a thermocycler (Electrocon corporation HBPX2220). 2.5 µl of the bisulfite treated DNA was used as template for the first round of PCR using external primers, then 2.5 µl of the PCR product was used as template for the nested PCR using internal primers. Primers were used at 0.2 µM and dNTP at 200 µM. Products were sequenced using the genoscreen sequencing facilities (http://www.genoscreen.fr/fr/). Chromatograms were edited using MEGA 6 [104] and DNA Chromatogram Explorer (http://www.dnabaser.com/download/chromatogram-explorer/). For each methylated cytosine, the percentage of methylation was determined according to five levels (0, 25, 50, 75, and 100%) depending on the height of the chromatogram peaks for the cytosine (non-converted methylated cytosine) compared to the height of the thymine (converted non-methylated cytosine) overlapping at the same position on the forward strand (and reversely on the reverse strand). For positive predictive value (PPV) calculation, the two regions: scaffold1255:365901.366597 and scaffold33832:13205.13768 were considered “diuron exposed profile” when individuals displayed percentage of methylation of 0 and above 25, respectively. The PPV value was calculated using a MedCalc based interface (https://www.medcalc.org/calc/diagnostic_test.php) and true positive took into account a positive occurrence in the two markers. The disease group was the 10 D samples; the non-disease group was the 10 SWC samples. True positives (a) are individuals for which DNA methylation profile is typical of “diuron exposed profile” in the D group. False positives (c) are individuals for whom DNA methylation profile is typical of “diuron exposed profile” but those individuals are from the control group. PPV equals a/(a + c) (Supplementary Fig. S3).

Expression Analysis of Target Genes

RNA extraction and cDNA synthesis were performed as previously described [40]. Briefly, individual whole tissue samples of oyster spat previously frozen in liquid nitrogen were homogenized in RL1 lysis buffer supplied by the NucleoSpin 96 RNA Tissue Core kit (Macherey-Nagel, France). Total RNA was extracted according to the manufacturer protocol. RNA samples were treated with DNase to prevent DNA contamination. cDNAs were synthesized from 300 ng of total individual RNA samples using VILO Superscript enzyme according to manufacturer’s instruction (Invitrogen). Primer specificity was tested performing a blastN against C. gigas genome version 9 (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Primer sequences and expected PCR product sizes are listed in Supplementary Table S1.

End point Polymerase Chain Reaction (PCR) was performed to check for the presence of alternative transcripts in the three candidates XLOC_005912 (CGI_10018082), XLOC_032353 (CGI_100 10582) and XLOC_032365 (CGI_10010585) and amplified, respectively, with the following primer pairs XLOC_005912.3for/XLOC_ 005912.5rev, XLOC_032353.2for/XLOC_032353.1rev and XLOC_0 32365.1for and XLOC_032365.3rev. A total of 20 individuals were tested per condition (SWC, SC and D). PCR amplification was performed in 25 µl using 0.625 U of GoTaq DNA (Promega kit: GoTaq G2 Hot Start Polymerase), 1× GoTaq buffer, 0.2 µM of each primer, 0.2 mM dNTPs, 2 mM MgCl2. PCR cycles were: 95 °C for 3 min, 35 cycles of 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 1 min and finally 72 °C for 10 min. PCR were accomplished on MastercycleEp Gradient Eppendorf thermocycler. 3 µl of the cDNA matrix were used. PCR product sizes were checked by electrophoresis (TAE 1%). For high resolution analysis on cDNA sequence, a 2.198 kb fragment from CGI_10018082 was amplified with XLOC _005912.3for and XLOC_005912.5rev. 1 µl of this PCR product was cloned into pCR4 (TOPO TA Cloning kit, Invitrogen) and the inserts (17 inserts per condition) derived from amplified cDNA of D (D3-17, D3-12, D3-19) and SWC (TE1-11, TE1-14, TE1-19) individuals were sequenced with primers (M13F, M13R and XLOC_005912.1rev) using the genoscreen sequencing facilities (http://www.genoscreen.fr/fr/). Chromatograms were edited using Sequencher Software, and a total of 34 sequences (combining several clones of each individual from each condition) were thoroughly examined for final analysis. Phylogeny analysis was performed using the on line phylogeny analysis tools [105].

Real time PCR (qPCR) were performed in duplicates in a total volume of 10 µl with 1.5 µl of H2O, 0.5 µl of each primers, 5 µl of MelTaq and 2.5 µl of cDNA (Roche kit: LightCycler® 480 SYBR Green I Master). qPCR were accomplished on the LightCycler480 II (Roche Diagnostics). For each reaction, the crossing point (Cp) was determined using the second derivative maximum method using Light Cycler Software version 3.3 (Roche Diagnostics). For each reaction, the level of transcription was normalized using the mean geometric transcription rate of three reference sequences: Elongation factor 1: EKC33063.1 Ubiquitin 60S ribosomal protein L40: EKC19428.1 40S ribosomal protein S6: AFJ91756.1 as previously described [106] (Supplementary Table S1). Difference between SWC and D individuals was tested statistically using Student’s test on R software.

Click here for additional data file.