Population epigenetic divergence exceeds genetic divergence in the Eastern oyster Crassostrea virginica in the Northern Gulf of Mexico

Abstract Populations may respond to environmental heterogeneity via evolutionary divergence or phenotypic plasticity. While evolutionary divergence occurs through DNA sequence differences among populations, plastic divergence among populations may be generated by changes in the epigenome. Here, we present the results of a genome‐wide comparison of DNA methylation patterns and genetic structure among four populations of Eastern oyster (Crassostrea virginica) in the northern Gulf of Mexico. We used a combination of restriction site‐associated DNA sequencing (RADseq) and reduced representation bisulfite sequencing (RRBS) to explore population structure, gene‐wide averages of F ST, and DNA methylation differences between oysters inhabiting four estuaries with unique salinity profiles. This approach identified significant population structure despite a moderately low F ST (0.02) across the freshwater boundary of the Mississippi river, a finding that may reflect recent efforts to restore oyster stock populations. Divergence between populations in CpG methylation was greater than for divergence in F ST, likely reflecting environmental effects on DNA methylation patterns. Assessment of CpG methylation patterns across all populations identified that only 26% of methylated DNA was intergenic; and, only 17% of all differentially methylated regions (DMRs) were within these same regions. DMRs within gene bodies between sites were associated with genes known to be involved in DNA damage repair, ion transport, and reproductive timing. Finally, when assessing the correlation between genomic variation and DNA methylation between these populations, we observed population‐specific DNA methylation profiles that were not directly associated with single nucleotide polymorphisms or broader gene‐body mean F ST trends. Our results suggest that C. virginica may use DNA methylation to generate environmentally responsive plastic phenotypes and that there is more divergence in methylation than divergence in allele frequencies.

a genome-wide comparison of DNA methylation patterns and genetic structure among four populations of Eastern oyster (Crassostrea virginica) in the northern Gulf of Mexico. We used a combination of restriction site-associated DNA sequencing (RADseq) and reduced representation bisulfite sequencing (RRBS) to explore population structure, gene-wide averages of F ST , and DNA methylation differences between oysters inhabiting four estuaries with unique salinity profiles. This approach identified significant population structure despite a moderately low F ST (0.02) across the freshwater boundary of the Mississippi river, a finding that may reflect recent efforts to restore oyster stock populations. Divergence between populations in CpG methylation was greater than for divergence in F ST , likely reflecting environmental effects on DNA methylation patterns. Assessment of CpG methylation patterns across all populations identified that only 26% of methylated DNA was intergenic; and, only 17% of all differentially methylated regions (DMRs) were within these same regions.
DMRs within gene bodies between sites were associated with genes known to be involved in DNA damage repair, ion transport, and reproductive timing. Finally, when assessing the correlation between genomic variation and DNA methylation between these populations, we observed population-specific DNA methylation profiles that were not directly associated with single nucleotide polymorphisms or broader genebody mean F ST trends. Our results suggest that C. virginica may use DNA methylation to generate environmentally responsive plastic phenotypes and that there is more divergence in methylation than divergence in allele frequencies.

K E Y W O R D S
Crassostrea virginica, DNA methylation, epigenetics, oyster, reduced representation bisulfite sequencing, salinity

| INTRODUC TI ON
In highly variable marine environments, phenotypic plasticity allows organisms to persist in suboptimum environmental conditions (Somero, 2010). On short time scales, this can be accomplished via environmentally responsive changes in gene expression (Evans & Hofmann, 2012). Over the medium to long term, there is increasing evidence that phenotypic plasticity may be modulated by changes in the epigenome (Eirin-Lopez & Putnam, 2018;Hofmann, 2017;Hu & Barrett, 2017;Verhoeven, vonHoldt, & Sork, 2016). These epigenomic changes allow for the persistent expression of a plastic phenotype that in some cases may be transmitted to subsequent generations (Rondon et al., 2017). As a result, there is significant interest in understanding how epigenetic variation correlates with genomic and environmental variation; and ultimately, how changes in the epigenome might facilitate the generation of environmentally responsive phenotypes. The push to address these questions has developed as a result of new sequencing techniques that provide high-resolution data on epigenomic traits and the necessary data for identifying when phenotypic effects of genotype-by-environment interactions are governed by changes in the epigenome (Trucchi et al., 2016;Van Gurp et al., 2016;Yaish, Peng, & Rothstein, 2014).
While there are multiple types of epigenetic mechanisms, 5mC DNA methylation (methylation at the 5th position in a cytosine molecule) has been the focus of a substantial amount of research to date using enzyme-linked immunoglobulin assays, methylation sensitive PCRs, and modified amplification polymorphism assays (Hu & Barrett, 2017;Kilvitis et al., 2014;Metzger & Schulte, 2016;Schrey et al., 2013). These methodologies are able to identify quantitative changes in DNA methylation; however, more targeted hypothesis testing is limited by a lack of high-resolution nucleotide-specific data among nonmodel marine invertebrates (Hofmann, 2017;Schrey et al., 2013). By contrast, reduced representation bisulfite sequencing (RRBS) methods provide nucleotide specific 5mC DNA methylation status using high-throughput short-read sequences on the Illumina platform (Trucchi et al., 2016;Van Gurp et al., 2016). In particular, the epiGBS (Van Gurp et al., 2016) method combines this RRBS approach with enzyme digestion that is not sensitive to methylation and therefore provides locus-specific methylation state across a large proportion of the genome, even when no reference genome is available.
Much of what is known about the functional role of epigenetic modifications in sessile marine invertebrates comes from research on the Pacific oyster (Crassostrea gigas). In C. gigas, DNA methylation appears to be an important molecular feature that directs developmental plasticity, gene regulation, and control over alternative splicing (Olson & Roberts, 2014;Riviere, 2014;Riviere et al., 2013;Song, Li, & Zhang, 2017). Specifically, increased levels of DNA methylation within intragenic regions in C. gigas are positively correlated with gene expression levels (Gavery & Roberts, 2013), while exon methylation is positively correlated with inclusion in mRNA transcripts (Song et al., 2017). Furthermore, changes in methylation patterns following parental exposure to environmental toxins appear to be heritable in the F1 generation (Rondon et al., 2017). Methylation research in C. gigas has been facilitated by the availability of a reference genome (Trucchi et al., 2016). The recent release of a reference genome for the Eastern oyster Crassostrea virginica (Gómez-Chiarri, 2018) provides an opportunity to explore similar patterns of DNA methylation with RRBS, and opens the door to future comparative epigenomic studies between these two closely related taxa.
The regulatory importance of epigenetic changes in the process of population divergence is increasingly being supported by models that suggest these nongenetic features influence the rate of adaptation to environmental stress (Klironomos, Berg, & Collins, 2016;Kronholm & Collins, 2016). Therefore, to better understand the genetic basis for population divergence, the interacting effects of genetic and epigenetic variation should be further explored.
Recent studies in wild populations have identified variation in DNA methylation that is associated with both acclimatization and broader population-specific epigenomic variation. These studies have identified methylation differences between freshwater and marine populations of the three-spined stickleback (Gasterosteus aculeatus) (Artemov et al., 2017) and between invasive and native populations of the pygmy mussel (Xenostrobus secures) (Ardura, Zaiko, Morán, Planes, & Garcia-Vazquez, 2017). Furthermore, multiple comparisons in species from clonal fish, tilapia, bats, and finches have found variation in DNA methylation that is greater than the standing genetic variation, lending support for the role of DNA methylation in alleviating the consequences of genotype-by-environment mismatches in a wide range of metazoans (Liu, Sun, Jiang, & Feng, 2015;Massicotte, Whitelaw, & Angers, 2011;Skinner et al., 2014;Wan et al., 2016).
Genotype-by-environment interactions shape population structure and can play an important role in establishing population specific DNA methylation patterns and in driving phenotypic differences between populations. Recent evidence has shown that reduced representation genomic sequencing methods such as restriction site-associated DNA sequencing (RADseq) are a valuable tool for evaluating population structure in marine invertebrates (Bernatchez et al., 2019;Silliman, 2019;Xuereb et al., 2018). The population differences that have been detected in many of these studies are surprising, as many marine invertebrates are broadcast spawners with large dispersal distances, a feature that should produce high gene flow across a species' range (Sanford & Kelly, 2011). However, previous assessments of population structure along the northern Gulf of Mexico (GOM) in Eastern oysters (Crassostrea virginica) have identified genetic variation between populations (Anderson, Karel, Mace, Bartram, & Hare, 2014). These genetic differences may be explained by either a shift in phenology between regions of the GOM or by the influence of genotype-by-environment mismatches removing some genotypes from any given estuary. In the northern GOM, one potential driver of genotype-environment mismatches is variation in salinity (Leonhardt, Casas, Supan, & Peyre, 2017).
In coastal Louisiana, populations of Crassostrea virginica are spread across estuaries that range from low salinity (<5 psu) to midhigh salinity (~20 psu) (Das et al., 2012). Within these estuaries, C. virginica plays a critical role as both a foundation species and an ecosystem engineer (Meyer, Townsend, & Thayer, 1997;Plunket & La Peyre, 2005). Previous research has shown that along the northern GOM, the interaction between temperature and salinity plays a major role in setting the species' environmental distribution (Eastern Oyster Biological Review Team, 2007;Shumway, 1996). While broad scale genotyping of oyster populations along the GOM has provided some evidence for genetic differentiation, there is little known about whether this genetic variation contributes to local adaptation to salinity (Varney, Galindo-Sánchez, Cruz, & Gaffney, 2009). In addition, regional salinity levels are also predicted to change dramatically in estuarine ecosystems due to changes in rainfall and anthropogenic alterations to hydrology (Das et al., 2012;Powell & Keim, 2015).
Previous research into the molecular responses of Crassostrea spp.
to environmental variation has identified significant responses to salinity at both the mRNA and protein levels (Chapman et al., 2011;Jones, Johnson, & Kelly, 2019). Therefore, examining the degree of genetic and epigenetic variation between populations of C. virginica that inhabit estuaries with distinct salinity conditions is an important step in understanding the observed signals of local adaptation and phenotypic plasticity in response to changes in environmental conditions along the northern GOM.
In this study, we characterize single nucleotide polymorphisms and DNA methylation patterns for four populations of Eastern oysters (Crassostrea virginica) along the northern GOM. This region is highly variable with respect to salinity, with shallow oyster reefs existing across regions with distinct annual salinity profiles ( Figure 1). Recent evidence suggests population-specific physiological responses to salinity between oysters from either high or low salinity reefs from this region (Leonhardt et al., 2017). In addition, we hypothesized that the Lake Fortuna population, located east of the Mississippi river, would be the most unique population both in terms of methylation and genomic variation. As such, we have used a combination of RADseq and epiGBS sequencing to further describe potential genetic and epigenetic drivers of these observed population-specific physiological responses.

| Collection sites
Oysters were collected from four sites that are known to differ in annual mean salinity. Salinities for each estuary were calculated as daily averages from June 1, 2006, to June 1, 2016, using salinity measurements (collected by either the United States Geological Survey or the Louisiana Department of Fish and Wildlife) from stations located as close to the oyster beds as possible ( Figure 1). These data were used to calculate the amount of time each reef spent at one of five salinity regimes known to impact oysters in the eastern GOM (which has higher average salinities than the northern GOM).

| Oyster collections
In order to characterize population-specific methylome patterns, 20 individual oysters from each site; Lake Calcasieu (LC), Lake Fortuna (LF), Sister Lake (SL), and Vermillion Bay (VB);were collected as part of a larger study between May 15, 2015, and February 10, 2016 ( Figure 1). Oyster from all sites was collected by dredging at a depth of ~3 m allowing for all oysters to have been collected from subtidal reefs. Salinity and temperature data were collected, and 10-year daily average is reported in addition to the average conditions for the 10 days preceding collection ( Figure 1c). Following collections, oysters were transported to Louisiana State University; gill tissue was dissected from each oyster, immersed in 95% ethanol, and stored at +4°C for 18-24 months. The remaining whole organisms were subsequently frozen and stored at −80°C. In order to gather morphological data, these oysters were thawed in September 2019. For individuals from LF and SL, individual specific length, wet weight, and presence/absence of gonads was assessed; however, due to freeze-thaw, it was not possible to identify the current sex of the individual. As oysters are sequential, protandric hermaphrodites, sexual stage has to be checked via gonads (Harding, Powell, Mann, & Southworth, 2013). For the other two sites (LC and VB), all 20 individuals were measured and weighed; however, the labeling of these individuals was lost, and so, only mean length and mean wet weight are reported.

| DNA extraction
DNA was extracted from all 80 individuals using the OMEGA E.Z.N.A. Tissue DNA Kit (D3396-01; Omega bio-tek) with a 2 min RNase A digestion to remove co-purified RNA. Extracted DNA purity was assessed on 260/280 and 260/230 ratios using a nanodrop spectrophotometer (ND1000; Thermo fisher Scientific). Presence of high molecular weight DNA was confirmed using a 1.5% agarose gel, and DNA concentration was verified using a Qubit 3.0 Flourometric dsDNA BR assay kit (Q32850; Life Technologies).
pl pipeline with the following settings "-p 2 --smooth --hwe -r 0.65 --min_maf 0.05 --bootstrap --bootstrap-reps 1000000 --structure --genepop --write_single_snp -vcf." These settings requires a locus to be represented by 65% of individuals in two of the four populations that have a minimum allele frequency of 0.05, F ST values are kernel-smoothed across chromosomes, and significance intervals based on resampling using 1,000,000 bootstraps (Appendix S1).
Population structure was investigated using the R program conStruct (v.1.0.3) and the population structure file generated from the Stacks program population. Samples were analyzed using the x.validation function in conStruct testing K from 1 to 4 using 90% of the dataset for training, 10 replicates, and 1,000 iterations per replicate (Bradburd, Coop, & Ralph, 2018). Input coordinates for each population are the same as those presented in Figure 1c.
The final K was chosen based on comparing layer contributions and predictive accuracy for models constructed both with and without geographic distances between each of the four estuaries ( Figure 2a,b).

| DNA methylation epiGBS library preparation
For epiGBS library preparation, a total of 500 ng of purified genomic
Trimmed, nondirectional, paired end reads were mapped to the reference genome, and CpG methylation was called using the software package bismark (v0.19.0) (Krueger & Andrews, 2011). The bismark commands used in the mapping allowed for 1 mismatch in a seed alignment of 10 with a minimum alignment score setting of −0.6 (--score_min L, 0, −0.6). These settings were selected to account for genomic variations between C. virginica collected from the northern GOM (this study) and the disease-resistant F I G U R E 2 (a) ConStruct crossvalidation analysis predictive accuracy of K from 1 to 4 using 90% of the dataset for training, 10 replicates, and 10,000 iterations per replicate. (b) ConStruct layer contributions of K from 1 to 4. (c) Population structure plot of K = 3 mapped to estuary locations. Downstream analysis of methylation focused on a K = 2 model that groups LC with SL and VB with LF inbred line from the U.S. East Coast used for the construction of the reference genome (Gómez-Chiarri, Warren, Guo, & Proestou, 2015). CpG methylation was extracted from the nondeduplicated mapped reads using the bismark command bismark_methylation_ extractor with the following commands; --ignore_r2 2, --bedGraph, --zero_based, --no_overlap, --cytosine_report, and -report.
Differential methylation was conducted on CpG features using the bismark coverage files and the R program MethylKit (v.1.2.4) (Akalin et al., 2012). Loci were first filtered using the filterByCoverage command to require a locus to have a low count of at least 10 reads to be included.
Methylation was then assessed in a pairwise bases using a 100 bp tiled approach with a step size of 100 bp and restricted to tiles that had been covered by at least eight individuals from each of the two populations. These parameters were chosen so that each 100-bp region would contain at least 1 CpG, would be approximately be equivalent to a single sequencing read (75 bp), and would be present in 40% of the population. Pairwise assessments of differential methylation were focused on contrasting the high salinity population with the other three populations. As such, the analysis described here reflects on both abiotic differences between estuaries over the 10-year daily mean differences between sites and the 10-day daily mean differences prior to sampling ( Figure 1c). Significantly differentially methylated regions

| Salinity variation among sites
Long-term salinity data for each site confirm variation in mean salinity, with Lake Calcasieu (LC) > Lake Fortuna (LF) > Sister lake (SL) > Vermillion Bay (VB; Figure 1). Examining distribution of time spent within each of the five salinity stress indices underscores the increasing frequency of time spent in low salinity conditions for VB and SL as compared with the two mid-high and high salinity sites ( Figure 1). Specifically, LC and LF have recorded salinities greater than 15 psu for over 50% of the year, while SL is below 15 psu for 71% of the year, and VB is below 3.5 for 57% of the year. These differences in exposure history support the use of these oyster reefs in our efforts to describe population-specific methylation profiles that may facilitate population differentiation and possibly local adaptation to low and high salinity in this region. however, this analysis only identified slightly more support for a K = 3 over K = 2 (Figure 2a,b). We interpret these results to suggest that while K = 3 is the best fit for the data, a K = 2 explains close to the same amount of population structure such that Lake Calcasieu and Sister Lake are highly similar and Lake Fortuna and Vermilion Bay are also highly similar populations (Figure 2).  (Figure 3).

| Differential methylation
The principle component analysis (PCA) revealed that the majority of samples between the four populations showed overlapping means ( Figure 4). Overall, 75.8% of all differentially methylated regions from all of the pairwise comparisons were located within gene bodies (promoters, introns, exons, 3′ and 5′ untranslated regions), 6.5% were within downstream regions, and 17.7% were located along intergenic regions (Figure 3). These distributions were tested against the randomly generated 100 bp regions and showed that the distribution of DMRs was higher than expected for introns, 3′ UTRs, and 5′ UTRs (p-value ≤ .012); lower than expected for exons and intergenic regions (p-value ≤ .008); but not significantly different from random for promoter (2 kb upstream) or downstream (2 kb downstream) regions (p-value > .05). Finally, gene ontology enrichment for each populationspecific differential methylation analysis was conducted for each pairwise comparison between DMRs that fell within the gene-body region using a Mann-Whitney U test with a background list of all genes for which there were sequence data within any of the comparisons.

F I G U R E 3 Distribution of percent methylation across genomic regions.
Global DNA methylation (black) and the mean distribution of randomly generated markers (gray) distribution of methylated cytosines across gene regions. Site specific mean methylation of DMRs for each genomic region is also shown in color; (Blue-Lake Calcasieu; Orange-Lake Fortuna; Green-Sister Lake; Purple-Vermillion Bay). Stars indicate samples are significantly different from the mean distribution of random markers

| Site-specific DNA methylation trends
Oysters from the Lake Calcasieu ( (Table 1). In addition, the most significant hypermethylated DMRs in the LC population were associated with genes coding for axonemal-like dynein heavy chain 8 gene, an atrial natriuretic peptide receptor 1-like isoform X3, and a serine/threonine-protein phosphatase subunit (Appendix S2). Gene ontology (GO) enrichment for DMRs that fell within gene bodies revealed enrichment in two biological process ontologies between LC and SL oysters.
These ontologies were carboxylic acid metabolic process and cellular amino acid metabolic process. GO enrichment also identified 2 molecular function ontologies between LC and VB which were identified as guanyl-nucleotide exchange factor activity and S-adenosylmethionine-dependent methyltransferase activity.
The Lake Fortuna site has the second highest annual mean salinity with approximately 51% of the year spent above 15 psu (Figure 1).
This site also differed from all others in timing of the sampling with individuals collected in May. As such, the short-term (10-day) difference in environmental conditions was largest for temperature  with SL. Hypomethylation of gene bodies in the VB population was associated with the gametocyte surface protein P230, a cytosolic carboxypeptidase 1-like protein, and an uncharacterized long noncoding RNA. Hypermethylated DMRs in the VB populations included genes associated with a germinal center associated nuclear protein-like, a cytoplasmic dynein 1 heavy chain 1-like isoform X1, and a palmitoyltransferase ZDHHC17-like isoform X1.
GO enrichment also identified significant enrichment between VB and both the LF and SL populations of the same broad category extracellular region. Genes associated with this ontology in both comparisons included regulatory genes such as WNT-5b, neurogenic locus notch homolog protein 4-like, and a phosphatidylinositol-glycan-specific phospholipase D-like gene (Appendix S2).

| Differential methylation between populations identified in conStruct
The final comparison of differences in methylation between populations used the results obtained from the conStruct analysis identified significant genomic structure at K = 3, with nearly equal support for K = 2. After assessment of the structure plot in Figure 2, we focused on exploring differences in methylation between a K = 2 model wherein the individuals sampled from Lake Calcasieu and Sister Lake were considered to be population 1 while the individuals sampled at Lake Fortuna and Vermilion Bay were considered to be population 2. This approach was chosen to test if the genomic population structure might explain the distribution of methylation between populations. Hierarchical clustering comparing population 1 with population 2 identified clustering based on population with only minimal overlap of individuals ( Figure 5a). Differences in the 10-day mean salinity and temperatures between the two population groups showed that population 1 was on average 4 psu higher than population 2, but population 2 was 7°C warmer than population 1. Global methylation patterns between these two meta-populations using a pairwise assessment of differential methylation identified 282 DMRs.
The top DMRs within gene bodies were hypomethylated in population 2 (LF + VB) when compared with population 1 (LC + SL) and were associated with genes coding for a histone acetyltransferase KAT6A-like gene, a anoctamin-7-like gene, a ATP-binding cassette sub-family A member 3-like gene, and a glutathione synthetaselike gene. The top hypermethylated DMRs were associated with a V-type proton ATPase subunit, a protogenin B-like gene, a sodium/hydrogen exchanger 9B2 gene, and dnaJ homolog subfamily C member 1-like gene ( Table 2). Gene ontology enrichment between these two populations revealed significant enrichment among hypermethylated genes associated with the ontology transferase activity, transferring glycosyl groups. Among the genes that were hypomethylated within the transferase activity was a glycylpeptide N-tetradecanoyltransferase 2-like gene that had a mean methylation of 7.7% methylation in population 2 and a mean methylation of 27.5% in population 1. This relationship was also true for the ATP-binding cassette sub-family A member 3-like gene and an uncharacterized protein.

| D ISCUSS I ON
There is an emerging interest in understanding the role DNA methylation plays in population responses to heterogeneous environments (Hofmann, 2017;Hu & Barrett, 2017;Metzger & Schulte, 2016). Recent physiological evidence suggests there is local adaptation between populations of C. virginica distributed across the northern Gulf of Mexico (GOM) (Leonhardt et al., 2017). Therefore, this study was designed to assess genetic and epigenetic variation between oysters collected from four coastal Louisiana sites spanning the range of salinity conditions inhabited by oysters in the GOM. By employing a combination of RADseq and comparative epigenomic analysis of oysters from these four sites, we were able to show significant population structure for genetic sequences, but much greater divergence in population specific methylation. This suggests that DNA methylation may play a greater role than sequence divergence in generating the previously observed phenotypic differences among C. virginica populations. While there are potential limitations to these data owing to the temporal distribution of collections and variation in prior exposure, these data provide important insight into epigenetic variation across known environmental gradients that may be contributing to the establishment of acclimatization and local adaptation among populations of C. virginica in this region.

| Genetic population structure
A surprising result from our analysis was the observed genetic similarities did not follow a geographic pattern, with the greatest similarities between Lake Fortuna and Vermillion Bay and between Lake Calcasieu and Sister Lake ( Figure 2). We had originally hypothesized that Lake Fortuna would be the most divergent population as the Mississippi river presents a significant barrier to larval dispersal along coastal Louisiana. Our findings may reflect the effects of human management.
In Louisiana, juvenile (seed) oysters are routinely moved from public seed grounds to private oyster leases. Public seed grounds have historically included Lake Calcasieu, Lake Fortuna, and Vermillion Bay, and there are private leases within or near Lake Fortuna, Vermillion Bay, and Sister Lake. In addition, hatchery produced larvae or juveniles have been outplanted at several of these locations by the Louisiana Department of Wildlife and Fisheries (Leonhardt, 2013). The influence of these oyster transplants and translocations may also explain the results from the conStruct analysis which suggests that the observed population structure is not driven by an isolation-by-distance effect.
The conStruct analysis also identified at least two distinct populations at each site, possibly reflecting the presence of both "native" and trans-

| Global DNA methylation patterns in C. virginica
The methylome sequencing approach used in this study revealed that approximately 14% of the C. virginica genome is methylated. This level of methylation is comparable to previous studies that report 15% of the C. gigas genome as being methylated (Gavery & Roberts, 2013;Olson & Roberts, 2014). The distribution of methylated regions across genomic features was also very similar to C. gigas, with a slightly higher percentage of regions within promoter regions being identified as methylated in C. virginica (Gavery & Roberts, 2013).
The observed overlap in frequency and distribution of methylation within gill tissue between these two closely related species suggests a conserved role for DNA methylation within oyster genomes. The range of population specific differentially methylated regions (38-282 regions) was consistently relegated to less than 10% of all measured regions (6.7%-8.6%), suggesting that the majority of methylation within the C. virginica genome is fixed and only a small portion is plastic. This is consistent with the concept that the major function of epigenomic modification is to control cellular differentiation and that environmentally induced de novo modifications are primarily important for maximizing phenotypic plasticity at both the cellular and organismal level (Eirin-Lopez & Putnam, 2018).
Therefore, focusing on the genomic regions that are associated with DMRs between the same tissue types from individuals collected at the same life-history stages will provide further evidence for which genes are important for maximizing genome-phenome diversity across a range of environmental conditions.
Among all pairwise assessment of differential methylation, gene ontology enrichment was only identified between comparisons of the low salinity site Vermillion Bay and both Lake Fortuna and Sister Lake.
In both cases, the enrichment was among the broad category extracellular region. Genes associated with this ontology were hypermethylated in Vermilion Bay and may reflect differences in developmental stages between individuals at these sites. One caveat regarding these data is that the collections presented here was measured on individuals that had significant variation in size; however, growth rate in C.virginica is dependent on environmental conditions; therefore, attributing age based on size is not a dependable method when comparing oysters from estuaries known to differ in salinity. However, presence or absence of gonads as was seen among individuals from Lake Fortuna suggests differences in reproductive status may have an influence on DNA methylation. Evidence for this can be seen in the pairwise assessments focused on the Lake Fortuna population that identified significant hypermethylation of genes associated with developmental processes.
Specifically, the hypermethylation of protogenin B and germinal center associated nuclear protein may reflect seasonal changes in DNA methylation that could be associated with reproductive timing. Future studies that account for these components will be crucial in validating the results presented here as samples varied in both time of collection and age of individuals.

| DNA methylation and population structure
Integrating the conStruct structure results from the RADseq analysis provided an important insight into the influence of genotype on DNA methylation. Specifically, adjusting the differential methylation comparisons to reflect a population structure of K = 2 identified 282 DMRs that showed enrichment for genes associated with transferring glycosyl groups. Specifically, enrichment was seen among hypomethylated genes coding for a gene with enzyme annotation of glycosyltransferase and the gene coding for a glutathione synthetase-like protein. These enzymes are important for the biosynthesis of many carbohydrates and detoxification that have been shown to be differentially regulated, at the transcript level, in response to stress in other Crassostrea species (Müller et al., 2018;Tanguy, Boutet, & Moraga, 2005). The hypermethylation of ketohexokinase, an enzyme involved in metabolism, suggests that differences in food availability between the estuaries may influence these methylation patterns. The hypermethylation of ion channels such as sodium/hydrogen exchanger 9B2 and the transient receptor potential cation channel suggests that the salinity differences between the two population groups may also influence differential methylation (Jones et al., 2019). Finally, the hypomethylated DMR located within the gene body of the ATP-binding cassette sub-family A member 3-like gene was found to be one of three regions that had less than 10% methylation in population 2. This gene is putatively involved in shell formation, and the transition from 7.5% methylation to 36.2% methylation may reflect physiologically significant changes in methylation (Feng, Li, Yu, Zhao, & Kong, 2015).
The relationship between genetic and epigenetic population structure suggests some influence of genetic structure on methylation state (see Figure 5). Despite this influence, the greater differences among populations in methylation suggests there is potentially both an environmental and developmental influence on epigenetic methylation status. While the RADseq analysis suggests that genetic structure may also be influenced by management related activities, the presence of a strong differential methylation of ion transporters and genes involved in stress responses provides evidence that populations have unique methylome profiles that may directly regulate the plasticity of these genes in response to variation in the environment.

| CON CLUS ION
Our analyses identified differential methylation among genes known to exhibit high plasticity in response to environmental conditions (sodium/hydrogen exchanger, glutathione synthetase), to developmental timing (protogenin B), and DNA damage repair (dnaJ homologs).
The variation in methylation around these genomic features potentially reflects both differences in timing of collection and differences in the environments the individuals were collected from. The DMRs identified here serve to highlight the diversity of methylation in the Eastern oyster and provide further evidence of the dynamic nature of CpG methylation in response to both environmental conditions and phenology. In addition, global methylation patterns revealed that relationships identified in C. gigas holds true for C. virginica and therefore supports the role of CpG methylation in phenotype environmental interactions for the genus Crassostrea. Assessing methylation patterns between populations identified using RADseq analysis also highlighted the effects of genomic structure on DNA methylation. Finally, our RADseq analysis also identified unexpected population structure potentially arising from management practices that may act counter to local adaptations to salinity within this coastal system. These results reinforce the need to combine genomic and epigenomic data when seeking to understand divergent population responses to heterogeneous environments.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data for this study will be made publicly available through the NCBI SRA database under accession PRJNA488288. STACKS code is available in Appendix S1, and all R code associated with this analysis is archived on Zenodo at https ://doi.org/10.5281/ zenodo.3551408.