A complete account of EDC exposures in the fish analyzed herein will be published elsewhere (DeCourten et al., unpublished). Briefly, embryos were obtained from an adult brood stock and processed according to methods in Porazinski et al. (2010), but without dechorionation. Fish [aged 8 hours post fertilization (hpf) – 21 days post hatch (dph)] in the parental (F0) generation were exposed to environmentally relevant concentrations of each of four chemicals separately: bifenthrin (Bif, 3.02 ng/L; Chem Services, West Chester, PA, United States; 99.5% pure mix of isomers), EE2 (6.79 ng/L; Sigma-Aldrich, St. Louis, MO, United States; CAS 57-63-6, >98% purity), levonorgestrel (Levo, 9.27 ng/L; USP, Rockville, MD, United States; 100% purity) and trenbolone (Tren, 9.60 ng/L; Spectrum Laboratory Products, Gardena, CA, United States; 100% purity). Concentrations for all chemicals are reported as average measured concentrations of exposure water, verified analytically using liquid chromatography-spectrometry described elsewhere (DeCourten et al., unpublished). Control groups were exposed to 10 μL/L of methanol to control for any effects of the solvent used for lipophilic chemicals and five experimental replicates were used for each treatment. New treatment water was mixed just prior to the daily 75% water changes. During periods in the 25 mL beakers (8 hpf-hatch at day 7–10) new water was mixed everyday by adding 10 μL (for BF or EE2) or 20 μL (for levonorgestrel or trenbolone) of EDC stock solutions in MeOH (0.1 mcg/mL). During exposure periods in the 1.4 L jars (hatch-21dph) water was mixed in 2–3 L batches by adding of 5 μL (bifenthrin and EE2) and 10 μL (levonorgestrel and trenbolone) for each liter made of EDC stock solutions in MeOH (1 mcg/L). Each replicate (n = 4–5) was maintained independently as described above throughout the course of the 3-generation study, with EDC exposures continuing until the 21-dph sampling time point (for the F0 generation only; Figure 1. After the initial 21-dph sampling time point of the F0 animals, a subset of animals was sampled for molecular endpoints, while the rest of the animals were transferred to clean water and no further EDC exposure occurred. At 21 dph remaining fish were transferred to 6 L glass containers where they were fed a diet of Hikari tropical fish food and live Artemia nauplii twice per day, with a 60% water change daily and reared to ∼120 dph. After 120 dph, fish were transferred to 20-gallon recirculating tanks where they were spawned in groups of ∼50 individuals (50:50 sex ratio) by placing strands of dye-free acrylic yarn into tanks overnight, allowing for spawning for roughly 3 h. After spawning, embryos were transferred to the laboratory where they were treated the same as F0 animals, but without chemical exposure, through the F2 generation. Thus, F1 animals were exposed to EDCs indirectly as primordial germ cells within F0 parents, while F2 animals were not exposed to EDCs at all. Any effects noted in the F1 animals would be considered “multigenerational” effects given that both F0 and F1 were either directly or indirectly exposed to the chemicals. Any effects noted in F2 animals, which were never exposed to the chemicals, are considered “transgenerational” effects.

DNA from fish in the EDC exposure experiment was extracted from two whole 21-dph larvae from each replicate tank (n = 4–5) using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). Samples were quantified using a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) and purity was assessed via electrophoresis on a 1% agarose gel and visualized on a Gel Doc^TM XR + Gel Documentation system (Bio-Rad, Hercules, CA, United States). After DNA extraction, approximately 500 ng of genomic DNA was treated with 1 μL of RNAseIF (NEB, Ipswich, MA, United States) for 30 min at 37°C followed by a 0.8X AmPure XP (Beckman Coulter, Carlsbad, CA, United States) clean up.

Genomic DNA was adjusted to a concentration of 1.0 ng/μl, and 1.25 ng of template gDNA was loaded on a Chromium Genome Chip. Whole genome sequencing libraries were prepared using Chromium Genome Library & Gel Bead Kit v.2 (10X Genomics, cat. 120258), Chromium Genome Chip Kit v.2 (10X Genomics, cat. 120257), Chromium i7 Multiplex Kit (10X Genomics, cat. 120262) and Chromium controller according to manufacturer’s instructions with one modification. Briefly, gDNA was combined with Master Mix, a library of Genome Gel Beads, and partitioning oil to create Gel Bead-in-Emulsions (GEMs) on a Chromium Genome Chip. The GEMs were isothermally amplified with primers containing an Illumina Read 1 sequencing primer, a unique 16-bp 10x bar-code and a 6-bp random primer sequence, and bar-coded DNA fragments were recovered for Illumina library construction. The amount and fragment size of post-GEM DNA was quantified prior using a Bioanalyzer 2100 with an Agilent High sensitivity DNA kit (Agilent, cat. 5067-4626). Prior to Illumina library construction, the GEM amplification product was sheared on an E220 Focused Ultrasonicator (Covaris, Woburn, MA, United States) to approximately 350 bp (55 s at peak power = 175, duty factor = 10, and cycle/burst = 200). Then, the sheared GEMs were converted to a sequencing library following the 10X standard operating procedure. The library was quantified by qPCR with a Kapa Library Quant kit (Kapa Biosystems-Roche) and sequenced on one lane of HiSeq4000 sequencer (Illumina, San Diego, CA, United States) with paired-end 150 bp reads. Using 420M sequencing reads from the genomic library, the Menidia genome was assembled using supernova (version 2.0.0) with default parameters, and annotated using GMAP with default parameters¹ and with the M. beryllina transcriptome (14,393 genes: Supplementary File “transcript_gene_list.txt”) established by Jeffries et al. (2015).

Reduced representation bisulfite sequencing (RRBS) libraries were generated using the Premium RRBS Kit from Diagenode (Liege, Belgium) following the manufacturer’s instructions. Fragment size distribution of resulting library pools was assessed via micro-capillary gel electrophoresis on a Bioanalyzer 2100 (Agilent, Santa Clara, CA, United States). The library pools were quantified by qPCR with a Kapa Library Quant kit (Kapa Biosystems/Roche, Basel Switzerland) and each pool was randomized and sequenced on six lanes of an Illumina HiSeq 4000 (Illumina, San Diego, CA, United States) run with single-end 100 bp reads.

Illumina reads were subjected to quality control using trim_galore² (version 0.4.5) under RRBS mode. Bases with quality higher than 30 were trimmed from the 3′ end of the reads first, followed by the removal of any adapter sequences from the reads. Reads less than 30 bp in length after trimming were discarded. Bismark (Krueger and Andrews, 2011) (version 0.19.0) with default parameters was used to map all reads that passed the quality control to the M. beryllina reference genome that was assembled using 10X linked read technology (not published). Methylated regions were extracted from the alignment using DMRfinder (version 0.3) (Gaspar and Hart, 2017) using default parameters (maximum 500 bp in length with no less than 3 CpG sites within). Differential methylation in each methylated region was analyzed using a binomial test. Multiple test correction was carried out using Benjamini–Hochberg procedure (Benjamini and Hochberg, 1995) which has been applied as an appropriate method for minimizing false positives in methylation datasets (Asomaning and Archer, 2012). We analyzed differential methylation by treatment (Bif, EE2, Levo, Tren) and generation (F0, F1, F2) relative to the control for that generation, and differential methylation was considered significant relative to the control if the Benjamini–Hochberg adjusted p-value was less than 0.05 (Supplementary Data Analysis). The genomic context of methylated regions with respect to genes of interest was generated by annotating the assembled genome using previously published transcriptomic data (Jeffries et al., 2015; Brander et al., 2016) using gmap (Wu and Watanabe, 2005). The relationship between a methylated region and a gene was generated by overlapping the methylated region with any defined genomic regions. We will focus our discussion on two different regions with respect to genes: (1) Upstream 1000 bp, which we considered generally representative of the gene promoter region and will be referred to as “promoter” for brevity, and (2) within gene body, which comprised the exonic gene region. While the true location of the promoter region varies by individual gene, others have also designated the promoter as 1000 bp upstream of the TSS (Brenet et al., 2011; Wan et al., 2016). Further, in a survey of eukaryotic genomes, vertebrate promoters were frequently identified in less than 1000 bp upstream of the corresponding gene (Yella et al., 2018).

To elucidate the effects of EDCs on differential methylation across generations, we performed three different types of analyses. First, we tracked differential methylation by treatment and generation at the level of the gene, to determine a snapshot of differentially methylated genes overall (promoter and gene body combined) and by region (promoter and gene body analyzed separately).

Next, we performed functional enrichment analyses to determine whether gene ontology (GO) terms or Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched in the treatment groups based on differential methylation. Genes of interest for each functional enrichment analysis were chosen based on a significant directional change of methylation (hyper, hypo) as well as gene region (promoter versus gene body) relative to controls. GO Term enrichment analyses for biological processes was performed using the R package topGO (version 2.37) (Alexa and Rahnenfuhrer, 2019). GO Term enrichment analyses for cellular component and molecular function were not included in the present work because we chose to focus on the biological processes that were affected by changes in methylation. KEGG Pathways Enrichment analyses were obtained using the R package KEGGREST (version 1.26.1) (Tenenbaum, 2020), based on the D. rerio pathways database from the R package org.Dr.eg.db (Carlson, 2019). D. rerio Ensembl ID’s were converted to KEGG IDs using the bioDBnet database³ (Mudunuri et al., 2009). Enriched GO terms or KEGG pathways were considered significant at the level of α = 0.01 (Fisher’s exact test) (Supplementary Data Analysis).

Finally, we focused on the differential methylation by gene region for a suite of 109 genes that have either been responsive to EDC exposure or those that have the potential to be responsive based on previous work and phenotypes associated with EDC disruption in M. beryllina (Jeffries et al., 2015; Frank et al., 2018; DeCourten et al., 2019a). These genes included hormone receptors, genes essential to osmoregulation and immune function, steroid metabolism, DNA methylation, oxidative stress (Supplementary Table S1).

The genome assembly included 52,786 scaffolds that sums up to 568,309,548 base pairs (bp) (Supplementary Figure S2). The N50 of the assembly was 4,365,106 bp, in 27 scaffolds. The completeness of the assembly in gene space was assessed using BUSCO (version 3.0.2) with the vertebrate database. Out of 2,586 core BUSCO genes, 82.9% were found to be complete, with 2,023 genes in single copy and 121 genes in duplicated copies. There were 224 fragmented BUSCO genes, and 218 BUSCO genes were missing from the assembly.

Analysis of RRBS data yielded the number of CpG sites between 1.7 million and 2.3 million in all samples. These CpG sites belonged to 315,410 methylated regions. The coverage information for each sample is summarized as a box plot (Supplementary Figure S3). The methylated regions that did not have more than 10 reads covered from at least 146 samples were removed before the differential methylation analysis. This filtering process reduced the number of methylated regions to 66,983.

All treatments (Bif, EE2, Levo, Tren) and generations (F0, F1, F2) showed evidence of differential methylation relative to the control in the analysis of all genes (14,393 genes informed by the M. beryllina transcriptome). The percentage of genes that were differentially methylated in their promoter and/or gene body ranged from 6% (Levo F0) to 11% (Bif F2) of all annotated genes (Figure 2). Differential methylation at the gene level was nominally higher in the F1and F2 generation animals (8–11% of all annotated genes) of each treatment compared with F0 animals (6–8% of all annotated genes) (Figure 2), suggesting a possible increase in multigenerational methylation effects (F1) and transgenerational methylation effects (F2) compared to direct parental effects (F0) caused by all four EDCs. A similar trend has been noted for other EDCs in transgenerational experiments. For an exposure occurring in the parental generation, oviparous fishes, including M. beryllina and zebrafish, the F2 generation would be the first generation to display transgenerational effect, compared the F3 generation in viviparous animals (Brander et al., 2017). In one study examining transgenerational epigenetic inheritance in gestating female rats exposed to the agricultural fungicide vinclozin, differential DNA methylation was compared between F3 and F1 generation sperm. Beck et al. (2017) found significantly more DMRs in F3 compared to F1 sperm, with no shared overlap between those DMRs, indicating distinct transgenerational epigenetic methylation patterns. The authors attributed the distinct methylation patterns between generations to the sensitive developmental period of direct exposure in F1s given that these animals were exposed during a period of PGC deprogramming and subsequent reprogramming. Another study on methylmercury-exposed zebrafish also showed more epimutations in F2 (transgenerational effect) compared with embryonically exposed F0 sperm, and most DMRs were unique between the two generations (Carvan et al., 2017). Thus, it seems that transgenerational epigenetic effects of EDCs can be highly variable, and dependent on epigenetic disruption during critical developmental stages. If epigenetic machinery (methylation, histones, ncRNAs, etc.) are disrupted during critical developmental windows in early generations, these perturbations appear to display an exacerbation of epigenetic changes in subsequent generations. Although the pattern of epigenetic erasure and reprogramming post fertilization and PGC differentiation that occurs in medaka (same phylogenetic group – Atherinomorpha – as Menidia) but not zebrafish (Jiang et al., 2013; Potok et al., 2013; Ortega-Recalde et al., 2019; Wang and Bhandari, 2019) has not been established in M. beryllina, evidence of transgenerational epigenetic effects as a result of dysregulation of epigenetic programming appears to be evident in both models (medaka, which shares the pattern with mammals and zebrafish, which does not) as well as in M. beryllina.

Aside from the evidence of epigenetic dysregulation driving differential methylation among larval M. beryllina, a subset of differential methylation effects (3–5% of all annotated genes, depending on treatment) were shared between F1 and F2 generations, sometimes originating in the F0 generation (Figure 2), and can be considered as candidates for evidence of transgenerational epigenetic inheritance in the traditional, strict sense. For example, Rissman and Adli (2014) discuss strict transgenerational epigenetic inheritance as comparable to imprinting, when a specific epimutation is established in the germ line and maintained in subsequent generations. The directional and positional methylation changes at the gene level among EDC treated M beryllina would need to be individually verified to confirm this phenomenon. Some genes (∼1% for each treatment) were differentially methylated in the F0 and F2 generations, but not in the F1, thus, the potential for those methylation changes to be directly inherited is unclear. Specific evidence for transgenerational epigenetic inheritance in the strict sense is discussed further below (see Differential Methylation in EDC Responsive Genes).

When distinguishing between gene regions (promoter, gene body), differential methylation was more frequent within the gene body compared with the promoter region for all treatments and generations (Figure 3). Despite our use of RRBS to detect differential methylation, which operationally favors promoter regions (Yong et al., 2016), we were still able to detect more gene body differential methylation than promoter methylation. CpG methylation in the promoter region has been correlated with decreased gene expression in fish (Wan et al., 2016). Methylation in the gene body, however, is not correlated with gene repression and has the potential to affect splicing (Laurent et al., 2010) and sometimes lead to increased gene expression (Jones, 2012). It is possible that our definition of the promoter region as only 1 kb upstream was too narrow to capture many of the CpG islands in functional promoters. Gene promoters vary in their proximity relative to TSS by individual gene, so any nominal cutoff for promoters will not accurately represent the functional promoter for any given gene. Further, gene body sizes are typically larger than promoters, which would favor gene body differential methylation detection in our study. Still, we captured a great amount of gene body methylation given our use of RRBS methods. Most differentially methylated genes showed clear evidence of either hypo- or hypermethylation relative to the control, but a minority (<6%) showed evidence of both hypo and hyper methylation in different loci within the same gene region (Figure 3), indicating that directional methylation changes were locally specific.

GO Term and KEGG pathway enrichment analyses use a standardized classification system and vocabulary to place genes into groups, making it easier to understand the overall functional implications of genes that, in this case, are differentially methylated more than would be expected from a random selection of genes. Our GO Term analysis focused on biological processes (and not molecular function or cellular component) while KEGG pathway analysis groups genes by their inclusion in a given pathway (Nguyen et al., 2019; The Gene Ontology Consortium, 2019). Functional analyses of GO term (biological process) enrichment based on differentially methylated gene regions yielded significant enrichment (α = 0.01) for a total of 144 GO terms when considering all treatments, generations, directional methylation, gene regions, while KEGG pathway analysis yielded 66 significantly enriched pathways (Supplementary Tables S4, S5). Differential gene body methylation produced a greater number of significantly enriched GO terms and KEGG pathways than did differential promoter methylation. As with overall differential methylation analyses, this skew may be related to the limited promoter window screened in the present study and the greater length in the gene body compared to the promoter (see above). Enrichment of the same term or pathway sometimes occurred across multiple treatments and generations (Table 1 and Figures 4, 5). It is important to note, however, that our use of a non-model organism with limited annotation for GO and KEGG enrichment analyses requires cautious interpretation, as enrichment results may be inflated by limited annotation information. Thus, we focus our discussion on enrichment of generalized categories that contain many enriched terms or those in which the number of differentially methylated annotated genes was high compared to the number of background genes expected in those gene ontology groups or pathways.

As part of our larger study, DeCourten et al. (unpublished) found craniofacial and/or skeletal deformities in at least one generation of all four EDC treatments. Our analysis of DNA methylation of enrichment of GO term and KEGG pathways produced changes in methylation for chondrocyte differentiation, cartilage condensation, embryonic morphogenesis, Wnt signaling (which is essential to embryonic development; Steinhart and Angers, 2018), and regulation of actin cytoskeleton (Figures 4, 5), overlapping with all treatments and generations for which craniofacial and/or skeletal deformity phenotypes were measured (Bif F0, EE2 F1, Levo F2, Tren F0 and F1) (DeCourten et al., unpublished). Wnt signaling is essential for embryonic development. Other studies have documented phenotypic growth and development effects caused by EDC exposure linked to changes in DNA methylation. For example, prenatal exposure of the EDC cadmium led to an overrepresentation of DNA methylation in genes that were essential for organ and morphological development and bone mineralization, although only in females (Kippler et al., 2013). Further, DNA methylation has been associated with craniofacial abnormalities (Alvizi et al., 2017). We also found altered cardiovascular deformities in Bif F0, EE2 F1, and Tren F0 animals, some of the same treatments for which we have identified corresponding DNA methylation enrichment in GO terms and KEGG pathways including cardiac conduction (TrenF0) and adrenergic signaling of cardiomyocytes (Bif F2, EE2 F1, Levo F1, and Tren F2 (Figures 4, 5), suggesting that altered methylation of cardiac genes and pathways may be linked to the physical deformities caused by EDC exposure. While we did not find perfect correspondence between the treatments/generations with cardiac deformities and the altered DNA methylation GO terms and KEGG pathways, the overlap suggests that DNA methylation may play a role in altered cardiac developmental phenotypes. Although information on the relationship between EDC exposure, cardiac phenotypes, and DNA methylation is lacking, early life exposure to bifenthrin and EE2 has been linked to altered cardiac phenotypes (Jin et al., 2009; Salla et al., 2016). The overlap between the physiological whole organism endpoints resulting from EDC exposure and the functional enrichment of related GO terms and KEGG pathways in our differential methylation analysis suggest that some of the multi- and transgenerational effects observed in DeCourten et al. (unpublished) may be the result of EDC altered epigenetic control, although further work would be needed to mechanistically confirm the connection.

We found evidence of a variety of methylation enriched GO terms and KEGG pathways that are critical to neurodevelopment (e.g., wnt-activated signaling pathway involved in forebrain neuron fate commitment, neural crest cell fate specification, embryonic neurocranium morphogenesis, mTOR signaling pathway, calcium signaling pathway; Figures 4, 5). One study (Frank et al., 2019) found that early life exposure to low (picomolar) concentrations of bifenthrin altered gene expression in calcium signaling pathways, and led to a latent olfactory predator avoidance cue behavioral phenotype in M. beryllina. The authors concluded that early-life effects of bifenthrin on neurodevelopment had effects later in life. Negative chemotaxis, or the movement away from chemicals, was the most frequently enriched GO term with differential methylation within the gene body across all treatments, affecting multiple generations within each treatment (Supplementary Table S1). It is possible that the alteration of methylation of genes involved in negative chemotaxis is yet another way in which EDCs have impacted neurodevelopmental processes. Others have documented other non-reproductive behavioral effects of estrogens and androgens in fish. Lagesson et al. (2019) exposed eastern mosquitofish (Gambusia holbrooki) to environmentally relevant levels (3 ng/L) of trenbolone at one of two temperature regimes for 21 days and found that trenbolone increased fish boldness behavior and altered predator escape behavior (but only at the higher temperature). Three-spined sticklebacks (Gasterosteus aculeatus) exposed to EE2 during early development and assayed for behavioral endpoints 8 months after exposures showed a reduction in anxiety behavior related to scototaxis (light/dark preference) (Porseryd et al., 2019). Zebrafish exposed to BPA (0.1 nM to 30 μM) as embryos showed evidence of hyperactivity in mid-range concentrations and altered transcription of genes involved in methylation (dnmt1 and cbs). The authors also found differential DNA methylation in many regions of the genome, primarily within gene bodies, and particularly in genes involved in neurodevelopment, suggesting the potential for differential DNA methylation of neurodevelopmental genes to affect observed differences in swimming behavior (Olsvik et al., 2019).

A number of GO terms and KEGG pathways involved in the regulation of the p53 tumor suppressor gene, generalized response to DNA damage, and regulation of the Wnt signaling pathway were enriched for at least one generation in each treatment. Indeed, Brander et al. (2016) showed that pathways involved in carcinogenesis were differentially expressed within M. beryllina exposed to bifenthrin. Further, trenbolone (Boettcher et al., 2011) and binary mixtures of tributyltin and EE2 (Micael et al., 2007) are known genotoxicants in fish. EDCs also play a role in dysregulating DNA methylation for genes important in carcinogenesis (Serman et al., 2014), a finding also supported elsewhere (Mirbahai et al., 2011). One study in mammals showed that most Wnt pathway gene bodies (and not promoters) were differentially methylated in cancer cells (Galamb et al., 2016). We found that Wnt pathway differential methylation was prominent in the gene bodies for all treatments, further strengthening the idea that EDCs may alter methylation in genes and pathways relevant to carcinogenesis. Gene body methylation of TP53 (which encodes p53) in somatic cells contributes to many types of cancer (Rideout et al., 1990; Jones, 2012). We found differential methylation in the gene bodies p53 pathways for all four EDCs. Exposure to EDCs has been widely linked with cancer phenotypes (Soto and Sonnenschein, 2010; Rachon, 2015; Karoutsou et al., 2017), and those cancer phenotypes linked to epigenetic changes (Bhandari, 2016), providing support that the differential methylation caused by the four EDCs in the present study may play a functional role in carcinogenesis.

Altered DNA methylation in some biological functions may explain the transgenerational transfer of altered methylation patterns. The enriched GO terms related to epigenetic control mechanisms (histone deubiquination, histone acetylation) occurred in F1 or F2 animals for three out of four EDC treatments, suggesting that genes important in epigenetic programming have been dysregulated in these treatments due to early life exposure during sensitive windows, thus potentially perpetuating EDC-induced epigenetic dysregulation in subsequent generations as in Beck et al. (2017) and Carvan et al. (2017).

Another mechanism potentially responsible for the differential DNA methylation that we observed in larval M. beryllina may be an EDC-induced change in metabolism, particularly glutathione metabolism. Glutathione conjugation is a common detoxification mechanism for estrogens (Zhu and Conney, 1998) and has been established as a primary pathway involved in bifenthrin (Brander et al., 2016), trenbolone (Evrard and Maghuin-Rogister, 1987), and levonorgestrel metabolism (Stanczyk and Roy, 1990). Adult M. beryllina exposed to 0.5, 5, and 50 ng L^–1 of bifenthrin showed a great deal of differential gene expression for metabolic processes (Brander et al., 2016), in line with our observed GO term and KEGG enrichment for differentially methylated metabolic processes and pathways. We found hypomethylated glutathione metabolism to be enriched in gene bodies for three of four EDC treatments (Bif, EE2, and Tren), which potentially affected glutathione metabolism. Evidence that metabolic alterations can contribute to epigenetic dysregulation is building. For example, more resources diverted to glutathione metabolism may compete with the resources needed by epigenetic machinery to methylate DNA and histones, thereby leading to a generalized hypomethylation (Lee et al., 2009; Oppold and Müller, 2017; Sharma and Rando, 2017). While we found both hyper- and hypomethylation as a result of EDC exposure, it is possible that some of the hypomethylation we observed was produced as a result of the demands of glutathione metabolism. Thus, the general glutathione pathway utilized to detoxify EDCs could be contributing to the differential methylation observed throughout the genome, and/or the differential methylation specific to the glutathione metabolism pathway may also be affecting the balance of methylation/demethylation in other parts of the genome.

The analysis of differential methylation in 109 potentially EDC-responsive genes (EDCRG) yielded 29 genes that were differentially methylated in at least one treatment/generation within their promoter and/or gene body regions. All treatments, generations, and gene regions (promoter, gene body) showed some differentially methylated EDCRGs, with the exception of the Bif F2 promoter region, for which no EDCRG differential methylation was noted (Supplementary Figures S6, S7). While a total of 18 EDCRG’s were differentially methylated in the promoter and/or gene body across all F0 treatments, 20 EDCRG showed a multigenerational (F1) differential methylation effect while 22 showed a transgenerational (F2) effect depending on treatment and generation (Supplementary Figures S6, S7).

Differentially methylated EDCRGs included those involved in calcium signaling (mtor, ryr1, calr, tgfb), osmoregulation (atnb233, atp1a1, atp1a1b, atp1a2), steroidogenesis (17ß-hsd, 3ß -hsd), hormone receptors (ar, esr2, esr3, gpr30), immune function and/or inflammation (cbr1, inha, stat1, tgfb), structure (krt13, pdlm1), DNA methylation (dnmt3a), metabolism (alas2, cyp1a1, smpd3), sex determination (dmrt1), protein translation (noa1), proteolysis and degradation (htra3, usp42), zinc transport (slc39a14), and more general signaling pathway regulation (ywhab) (Supplementary Figures S6, S7). By frequency, the genes with the most alterations to their methylation patterns across all treatments and generations were dmrt1, ar, 17ß-hsd, atp1a1, atnb233, atp1a1b, ywhab, gpr30, cbr1, mtor, atp1a2, htra3 (Supplementary Table S8). The differential methylation we observed in these genes showed concordance with the generalized patterns of altered methylation enrichment in GO term and KEGG pathway analyses. Many of those genes showed differential methylation across generations within the same treatment, primarily within the gene body. Differential methylation occurring across generations was sometimes consistent (e.g., hypomethylation in the promoter of dmrt1 for Tren F0 through F2 generations), although directional methylation often alternated by generation (e.g., hypomethylation of atp1a1 within the gene body in EE2 F0 and F2, but hypermethylation of the same gene in F1), with hypomethylation in the gene body prominent in generations F0 and F2, but hypermethylation often prominent in the F1 generation (Figure 6). As with the overall analysis of all genes, a subset of EDCRGs showed patterns of both hyper- and hypomethylation, within the same greater region (promoter, gene body), although a singular directional differential methylation was the most common.

The majority of differential gene methylation documented for EDCRGs did not occur at the same sites within a given gene across the generations. Or, if they did, the direction of methylation change relative to the control often alternated (Supplementary Table S8), potentially indicating the interplay of epigenetic feedback loops (e.g., Bonasio et al., 2010). However, the same directional methylation changes (hypo, hyper) at the same loci were sometimes noted in F0 and F1 generations of the same treatment (17ß-hsd in Tren; ar in Bif, EE2, and Levo; atp1a1 in Bif; cbr1 in Bif, EE2, and Levo; dmrt1 in Bif and EE2) indicating identical effects occurring from exposure in both generations or potentially a multigenerational inheritance pattern. Although this phenomenon was less prevalent in F1 and F2 generations of the same treatment (17ß-hsd in Bif; dmrt1 in Tren; htra3 in Bif and Levo), the maintenance of specific methylation from an exposed generation (F1) to an unexposed generation (F2) indicates a strict form of transgenerational epigenetic inheritance (Figure 6). These examples of transgenerational epigenetic inheritance provide evidence of germline-established differential methylation, in contrast to the evidence of epigenetic programming dysregulation that we also found.

Doublesex And Mab-3 Related Transcription Factor 1 (dmrt1) was one of the genes that showed evidence of transgenerational epigenetic inheritance and was the most frequently differentially methylated gene across all treatments and generations. Dmrt1 is an important gene in sexual determination across many species, including those that exhibit temperature dependent sex determination, as Menidia species do (DeCourten et al., 2017; Huang et al., 2017) and in the closely related Japanese medaka (Otake et al., 2008). In the marine half-smooth tongue sole (Cynoglossus semilaevis), the level of dmrt1 demethylation was responsible for male gonad development (Shao et al., 2014). It is possible that the skewed sex ratios that are common in EDC exposed populations (Jobling et al., 1998; Kidd et al., 2007; Hua et al., 2015; Orn et al., 2016; DeCourten et al., 2017) may be the result of epigenetic reprogramming at the dmrt1 locus, as evidenced from the differential methylation caused by both estrogenic and androgenic EDCs in the present study.

We found differential methylation in genes involved in calcium signaling (calr, ryr1, mtor, tgfb) depending on treatment and generation, although all four EDCs yielded at least one calcium signaling gene with differential methylation (Supplementary Figures S6, S7). Alterations in calcium signaling gene expression as a result of bifenthrin exposure in M. beryllina and zebrafish have been established (Frank et al., 2018, 2019), and in response to other EDCs as well (Brander, 2013). Early life exposures of EE2 in coho salmon (Oncorhynchus kisutch) led to changes in transcription for both tgfb and calcium signaling pathways (Harding et al., 2013). The functional significance of the altered methylation in calcium signaling genes requires further study, but our preliminary results suggest that DNA methylation may play a regulatory role in altered calcium signaling pathway function in EDC-exposed fish.

We found differential methylation in the gene body of dnmt3a, the DNA methyl transferase responsible for de novo methylation (Okano et al., 1999; Hermann et al., 2004), in three different treatments, thus providing a potential mechanism to help explain DNA methylation dysregulation in exposed animals, although more work would need to be performed to confirm the functional relevance of this differential methylation. Other EDCs (TCDD, DES, PCB153) have been shown to alter transcription of Dnmts, which may in turn have direct effects on methylation (Wu et al., 2006).

Hormone receptors and regulators of steroidogenesis play an important role in development and reproduction. A recent multigenerational study with M. beryllina showed that exposure to 1 ng/L EE2 produced gene expression (gpr30, 17ß-hsd) changes in directly exposed generations of fish. Exposure to 1 ng/L bifenthrin instead produced latent effects in indirectly exposed F1 generation fish (GPR30) (DeCourten et al., 2019a). We found evidence of differential methylation for gpr30 in EE2 F0 and F2 and Bif F1 animals, while 17ß-hsd was differentially methylated in nearly all treatments and generations, indicating that a change in methylation state for these genes may be related to gene expression differences observed previously (DeCourten et al., 2019a). We also documented frequent differential methylation of the ar, another hormone receptor, across most treatments and generations. Casati et al. (2013) found that ar expression was modulated by other EDCs (PCBs), and that the altered expression had a basis in epigenetic mechanisms. Changes in 17ß-hsd, gpr30, and ar methylation have been correlated with cancer (Bhavani et al., 2008; Tian et al., 2012; Manjegowda et al., 2017), thus there is the potential for differential methylation of these genes to have significant phenotypic effects during both development and later in life.

A limited analysis of differential gene methylation in larval M. beryllina from the same treatments as the present study showed little correlation between changes in gene expression and DNA methylation in a suite of twenty genes (measured at 21 dph). In fact, only one gene, 17ß-hsd, showed hypomethylation within the gene body and significantly decreased expression in EE2 F0 animals (DeCourten et al., unpublished). The lack of concordance between gene expression changes and methylation could indicate that the differential methylation we observed is not functionally relevant and does not, therefore, translate to altered gene expression. Or, the lack of correlation could be a result of the single time point at which gene expression was measured, with epigenetic modifications being more stable than gene expression responses to EDCs. It is also possible that our approach of tracking methylation in the whole body larval fish, rather than specific tissues or cell-types created a dilution effect, making it difficult to track and correlate gene expression and DNA methylation. Future studies involving specific cell types (i.e., gonad) may allow for more prominent differential methylation trends to be revealed on a more mechanistically relevant basis. Other studies have also documented a lack of correspondence between DNA methylation and gene expression, even with tracking tissue-specific methylation (Aluru et al., 2018; Ryu et al., 2018). Further work is needed to determine the functional relationship between EDC-altered gene expression and differential methylation, ideally with multiple timepoint measurements and tissue and/or cell-specific DNA methylation profiling. Still, our detection of differential DNA methylation in all treatments and generations relative to the control provides evidence that EDCs drive changes in methylation on a multi- and transgenerational scale – a phenomenon that has already been established in mouse and human models (Susiarjo et al., 2007).