Epigenetic modification, specifically DNA methylation, is one possible mechanism for intergenerational plasticity. Before inheritance of methylation patterns can be characterized, we need a better understanding of how environmental change modifies the parental epigenome. To examine the influence of experimental ocean acidification on eastern oyster (Crassostrea virginica) gonad tissue, oysters were cultured in the laboratory under control (491 ± 49 μatm) or high (2550 ± 211 μatm) pCO₂ conditions for 4 weeks. DNA from reproductive tissue was isolated from five oysters per treatment, then subjected to bisulfite treatment and DNA sequencing. Irrespective of treatment, DNA methylation was primarily found in gene bodies with approximately 22% of CpGs (2.7% of total cytosines) in the C. virginica genome predicted to be methylated. In response to elevated pCO₂, we found 598 differentially methylated loci primarily overlapping with gene bodies. A majority of differentially methylated loci were in exons (61.5%) with less intron overlap (31.9%). While there was no evidence of a significant tendency for the genes with differentially methylated loci to be associated with distinct biological processes, the concentration of these loci in gene bodies, including genes involved in protein ubiquitination and biomineralization, suggests DNA methylation may be important for transcriptional control in response to ocean acidification. Changes in gonad methylation also indicate potential for these methylation patterns to be inherited by offspring. Understanding how experimental ocean acidification conditions modify the oyster epigenome, and if these modifications are inherited, allows for a better understanding of how ecosystems will respond to environmental change.

Introduction

As increased anthropogenic carbon dioxide is expected to create adverse conditions for calcifying organisms (IPCC, 2019), efforts have been made to understand how ocean acidification impacts ecologically and economically important organisms like bivalves (Parker et al., 2013; Ekstrom et al., 2015). Bivalve species are sensitive to reduced aragonite saturation associated with ocean acidification, with larvae being particularly vulnerable (Barton et al., 2012; Waldbusser et al., 2014). Shell structure may be compromised in larvae, juveniles, and adults (Gazeau et al., 2007; Kurihara et al., 2007; Beniash et al., 2010; Ries, 2011). Aside from affecting calcification and shell growth, ocean acidification can impact protein synthesis, energy production, metabolism, antioxidant responses, and reproduction (Tomanek et al., 2011; Timmins-Schiffman et al., 2014; Dineshram et al., 2016; Boulais et al., 2017; Omoregie et al., 2019).

Additionally, adult exposure to ocean acidification may impact their larvae [reviewed in Ross et al. (2016), Byrne et al. (2019)]. For example, adult Manila clams (Ruditapes philippinarum) and mussels (Musculista senhousia) reproductively conditioned in high pCO₂ waters yield offspring that exhibit significantly faster development or lower oxidative stress protein activity in those same conditions (Zhao et al., 2018, 2019). In contrast, northern quahog (hard clam; Mercenaria mercenaria) and bay scallop (Argopecten irradians) larvae may be more vulnerable to ocean acidification and additional stressors when parents are reproductively conditioned in high pCO₂ waters (Griffith and Gobler, 2017). Some species exhibit both positive and negative carryover effects [e.g., Saccostrea glomerata; (Parker et al., 2012, 2017)]. Intergenerational effects have also been documented when adult exposure to ocean acidification does not coincide with reproductive maturity [e.g., Crassostrea gigas; (Venkataraman et al., 2019)]. Although intergenerational carryover effects are now at the forefront of ocean acidification research in bivalve species, the mechanisms responsible for these effects are still unclear.

Epigenetics is the next frontier for understanding how environmental memory may modulate phenotypic plasticity across generations (Eirin-Lopez and Putnam, 2018). Epigenetics refers to changes in gene expression that do not arise from changes in the DNA sequence, with methylation of cytosine bases being the most studied mechanism (Bird, 2002; Deans and Maggert, 2015). Unlike highly methylated vertebrate genomes, marine invertebrate taxa have sparse methylation throughout their genomes, similar to a mosaic pattern (Suzuki and Bird, 2008). Genes that benefit from stable transcription, such as housekeeping genes, tend to be more methylated, while environmental response genes that are less methylated are prone to more spurious transcription and alternative splicing patterns, thereby possibly increasing phenotypic plasticity (Roberts and Gavery, 2012; Dimond and Roberts, 2016; Gatzmann et al., 2018). Increased levels of DNA methylation can also correlate with increased transcription. Several base pair resolution studies in C. gigas demonstrate a positive association between DNA methylation and gene expression that is consistent across cell types (Roberts and Gavery, 2012; Gavery and Roberts, 2013; Olson and Roberts, 2014). Since DNA methylation could provide a direct link between environmental conditions and phenotypic plasticity via influencing gene activity, elucidating how invertebrate methylomes respond to abiotic factors is crucial for understanding potential acclimatization mechanisms (Bossdorf et al., 2008; Hofmann, 2017).

While bivalve species have been used as model organisms to characterize marine invertebrate methylomes, how ocean acidification affects bivalve DNA methylation is poorly understood. Methylation responses to ocean acidification have been studied in multiple coral species. When placed in low pH conditions (7.6–7.35), Montipora capitata did not demonstrate any differences in calcification, metabolic profiles, or DNA methylation in comparison to clonal fragments in ambient pH (7.9–7.65) (Putnam et al., 2016). DNA methylation increased in another coral species, Pocillopora damicornis, in addition to reduced calficiation and more differences in metabolic profiles (Putnam et al., 2016). The coral Stylophora pistillata also demonstrates increased global methylation as pH decreases (pH_treatment = 7.2, 7.4, 7.8; pH_control = 8.0), with methylation reducing spurious transcription (Liew et al., 2018b). Combined whole genome bisulfite sequencing and RNA sequencing revealed differential methylation and expression of growth and stress response pathways controlled differences in cell and polyp size between treatments (Liew et al., 2018b). The association between DNA methylation and phenotypic differences in these corals demonstrates that epigenetic regulation of genes is potentially important for acclimatization and adaptation to environmental perturbation. Recent examination of C. virginica methylation patterns in response to a natural salinity gradient suggests that differential methylation may modulate environmental response in this species (Johnson and Kelly, 2019).

There is evidence that suggests that methylation patterns can be inherited in marine invertebrates. For example, purple sea urchin (Strongylocentrotus purpuratus) offspring have methylomes that reflect maternal rearing conditions (Strader et al., 2019). Different parental temperature and salinity regimes influence larval methylomes in Platygyra daedalea (Liew et al., 2018a). In the Pacific oyster (C. gigas), parental exposure to pesticides influence DNA methylation in spat, even though the spat were not exposed to these conditions (Rondon et al., 2017). Methylation changes in gametes are likely the ones that could be inherited, and may play a role in carryover effects. Before determining if DNA methylation is a viable mechanism for altering the phenotypes of offspring or subsequent generations, the epigenome of bivalve reproductive tissue in response to ocean acidification must be characterized.

The present study is the first to determine if ocean acidification induces differential methylation in reproductive tissue in the eastern oyster (Crassostrea virginica). Adult C. virginica were exposed to control or elevated pCO₂ conditions. We hypothesize that ocean acidification will induce differential methylation in C. virginica gonad tissue, and that genes with differentially methylated loci will have biological functions that could allow for acclimatization to environmental perturbation. Understanding how experimental ocean acidification conditions modify the oyster epigenome, and if these modifications are inherited, allows for a better understanding of how ecosystems will respond to environmental change.

Materials and Methods

Experimental Design

Adult C. virginica (9.55 cm ± 0.45) were collected from an intertidal oyster reef in Plum Island Sound, MA, United States (42.681764, −70.813498) in mid-July 2016. The oysters were transported to the Marine Science Center at Northeastern University (Nahant, MA, United States), where they were cleaned and randomly assigned to one of six flow-through tanks (50 L) maintained at ambient seawater conditions. Oysters were acclimated for 14 days under control conditions (500 μatm; 14–15°C) before initiating a 28-day experimental exposure. Half of the tanks remained at control pCO₂ conditions (500 μatm, Ω_calcite > 1), while the other half were ramped up to elevated pCO₂ conditions (2500 μatm, Ω_calcite < 1) over 24 h. This elevated treatment is consistent with observations in other estuarine ecosystems that oysters inhabit (Feely et al., 2010), although pH in nature only stays as extreme for short periods of time (e.g., hours). Moreover, the extreme treatment was also chosen to increase precision and therefore power to detect a response (Whitlock and Schluter, 2014).

Treatment conditions were replicated across three tanks, with oysters distributed evenly among tanks (1–2 oysters per tank). Each tank had an independent flow-regulator that delivered fresh, natural seawater at approximately 150 ml min^–1. Carbonate chemistry was maintained independently for each tank by mixtures of compressed CO₂ and compressed air at flow rates proportional to the target pCO₂ conditions. Gas flow rates were maintained with Aalborg digital solenoid-valve-controlled mass flow controllers (Model GFC17, precision = 0.1 mL/min). Within a treatment, tanks were replenished with fresh seawater and each tank was independently bubbled with its own mixed gas stream, with partial recirculation and filtration with other tanks in the treatment. As a result, the carbonate chemistry (i.e., the independent variable by which the treatments were differentiated) of the replicate tanks were slightly different from each other, which is evidence of their technical independence. Temperature was maintained at 15°C using Aqua Euro United States model MC-1/4HP chillers coupled with 50-watt electric heaters. Average salinity was determined by the incoming natural seawater and reflected ambient ocean salinity of Massachusetts Bay near the Marine Science Center (Latitude = 42.416100, Longitude = −70.907737).

Oysters were fed 2.81 mL/day of a 10% Shellfish Diet 1800 twice daily following Food and Agriculture Organization’s best practices for oysters (Helm and Bourne, 2004). Five oysters were collected from each treatment at the end of the 28-day exposure. They were immediately dissected with gonadal tissue harvested and immediately flash frozen. Partial gamete maturation was evident upon visual inspection.

Measurement and Control of Seawater Carbonate Chemistry

The carbonate chemistry of tanks was controlled by bubbling mixtures of compressed CO₂ and compressed air at flow rates proportional to the target pCO₂ conditions. The control pCO₂ treatments were maintained by bubbling compressed ambient air only.

Temperature, pH, and salinity of all replicate tanks was measured three times per week for the duration of the experiment. Temperature was measured using a glass thermometer to 0.1°C accuracy, pH was measured using an Accumet solid state pH electrode (precision = 1 mV), salinity was measured using a YSI 3200 conductivity probe (precision = 0.1 ppt). Every 2 weeks, seawater samples were collected from each replicate tank for analysis of dissolved inorganic carbon (DIC) and total alkalinity (A_T). Samples were collected in 250 mL borosilicate glass bottles sealed with a greased stopper, immediately poisoned with 100 μL saturated HgCl₂ solution, and then refrigerated. Samples were analyzed for DIC via coulometry and Alk_T via closed-cell potentiometric Gran Titration with a VINDTA 3C (Marianda Corporation). Other carbonate system parameters, including Ω_calcite, pH, and pCO₂, were calculated from DIC, A_T, salinity, and temperature using CO2SYS software version 2.1 (Lewis and Wallace, 1998; Van Heuven et al., 2011), using the seawater pH scale (mol/kg-SW) with K1 and K2 values from Roy et al. (1993), a KHSO₄ value from Dickson (1990), and a [B]_T value from Lee et al. (2010).

MBD-BS Library Preparation

DNA was isolated from five gonad tissue samples per treatment using the E.Z.N.A. Mollusc Kit (Omega) according to the manufacturer’s instructions. Isolated DNA was quantified using a Qubit dsDNA BR Kit (Invitrogen). DNA samples, ranging from 12.8 to 157 ng/μL, were placed in 1.5 mL centrifuge tubes and sonicated using a QSONICA CD0004054245 (Newtown, CT) in 30 s interval periods over 10 min at 4°C and 25% intensity. Shearing size (350 bp) was verified using a 2200 TapeStation System (Agilent Technologies). Samples were enriched for methylated DNA with the MethylMiner kit (Invitrogen). A single-fraction elution using 400 μL of high salt buffer was used to obtain captured DNA. After ethanol precipitation, 25 μL of buffer was used for the final elution. Library preparation and sequencing was performed by ZymoResearch using Pico Methyl-Seq Library Prep Kit (Cat. #D5455). Libraries were then barcoded and pooled into two lanes (eight samples in one and two in another) to generate 100bp paired-end reads on the HiSeq1500 sequencer (Illumina, Inc.).

Global Methylation Characterization

Sequences were trimmed with 10 bp removed from both the 5′ and 3′ ends using TrimGalore! v.0.4.5 (Martin, 2011). Quality of sequences was assessed with FastQC v.0.11.7 (Andrews, 2010). The C. virginica genome (NCBI Accession GCA_002022765.4) was prepared using Bowtie 2-2.3.4 [Linux x84_64 version; (Langmead and Salzberg, 2012)] within the bismark_genome_preparation function in Bismark v.0.19.0 (Krueger and Andrews, 2011). Trimmed sample sequences were then aligned to the genome using Bismark v.0.19.0 (Krueger and Andrews, 2011) with non-directionality specified and alignment score set using -score_min L,0,-1.2. Alignment files (i.e., bam) were deduplicated (deduplicate_bismark), sorted and indexed using SAMtools v.1.9 (Li et al., 2009). Methylation calls were extracted from deduplicated files using bismark_methylation_extractor.

Various C. virginica genome feature tracks were created for downstream analyses using BEDtools v2.26.0 (Quinlan and Hall, 2010). Genes, mRNA, coding sequences, and exons were derived directly from the C. virginica genome on NCBI (Gómez-Chiarri et al., 2015). The complement of the exon track was used to identify introns, and coding sequences were subtracted from exons to identify untranslated regions of exons (UTR). Exon locations were removed from the complement of the gene track to define intergenic regions. Putative promoter regions were defined as the 1 kb upstream of transcription start sites. Putative transposable elements were identified using RepeatMasker (v4.07) with RepBase-20170127 and RMBlast 2.6.0 (Smit et al., 2013; Bao et al., 2015). All species available in RepBase-20170127 were used to identify transposable elements.

Overall C. virginica gonad methylation patterns were characterized using information from all samples. Individual CpG dinucleotides with at least 5× coverage in each sample were classified as methylated (≥50% methylation), sparsely methylated (10–50% methylation), or unmethylated (<10% methylation). The locations of all methylated CpGs were characterized in relation to putative promoter regions, UTR, exons, introns, transposable elements, and intergenic regions. We tested the null hypothesis that there was no association between the genomic location of CpG loci and methylation status (all CpGs vs.methylated CpGs) with a chi-squared contingency test (chisq.test in R Version 3.5.0).

Methylation islands were determined to characterize overall methylation in the C. virginica genome using a sliding window analysis based on (Jeong et al., 2018). Islands were defined as areas of the genome with enriched levels of methylated CpGs (>50% methylation). To define methylation islands, each chromosome was examined using an initial 500 bp window starting at the first methylated CpG. If the proportion of methylated CpGs in the window was greater than 0.2, the window was extended by 50 bp; if not, the analysis proceeded to the next methylated CpG. Windows were continually extended until the proportion of methylated CpGs in the window fell below the 0.2 criteria. The location of methylation islands in the genome were characterized using BEDtools intersect v2.26.0.

Differential Methylation Analysis

Differential methylation analysis for individual CpG dinucleotides was performed using methylKit v.1.7.9 in R (Akalin et al., 2012) using deduplicated, sorted bam files as input. Only CpGs with at least 5× coverage in each sample were considered for analysis. Methylation differences between treatments were obtained for all loci in the CpG background using calculateDiffMeth, a logistic regression built into methylKit. The logistic regression models the log odds ratio based on the proportion of methylation at each locus:

$$log\left(\frac{{\mathrm{\pi }}_i}{1-{\mathrm{\pi }}_i}\right)={\mathrm{\beta }}_0+{\mathrm{\beta }}_{1}Treatmen{t}_i$$

A differentially methylated locus (DML) was defined as an individual CpG dinucleotide with at least a 50% methylation change between treatment and control groups, and a q-value < 0.01 based on correction for false discovery rate with the SLIM method (Wang et al., 2011). Hypermethylated DML were defined as those with significantly higher percent methylation in oysters exposed to high pCO₂ conditions, and hypomethylated DML with significantly lower percent methylation in the high pCO₂ treatment. A Principal Components Analysis (PCA) was performed for differentially methylated loci (DML) for oyster sample methylation profiles between treatments, then compared to a PCA for all MBD-enriched CpG loci. The location of DML were characterized in relation to putative promoter regions, UTR, exons, introns, transposable elements, and intergenic regions using BEDtools intersect v2.26.0. Loci that did not overlap with the aforementioned genomic features were also identified. A chi-squared contingency test was used to test the null hypothesis of no association between genomic location and methylation status between MBD-enriched CpGs and DML. To describe the location of DML across different gene architectures, the position of DML in the gene was scaled from 0 to 100 bp.

Enrichment Analysis

Functional enrichment analyses were used to determine if any biological processes were overrepresented in genes based on individual CpG methylation levels. Enrichment was conducted with GO-MWU, a rank-based gene enrichment method initially developed for analyzing transcriptomics data (Wright et al., 2015). Instead of only using genes with DML, GO-MWU identifies GO categories that are overrepresented by genes with any CpGs, allowing for more data to contribute to any trends. GO-MWU scripts and a gene ontology database were downloaded from the GO-MWU Github repository¹.

A gene list and table of significance measures were used as GO-MWU analysis inputs. The gene list contained Genbank IDs and all associated gene ontology terms. For the table of significance measures, Genbank IDs were matched with the smallest P-value for associated CpGs analyzed by methylKit. To match the Genbank IDs to CpG loci within mRNAs and create the gene list, overlaps between the C. virginica mRNA track from NCBI and the CpG background used in methylKit were obtained using BEDtools intersect v2.26.0. The mRNAs were then annotated with Uniprot Accession codes using a BLASTx search [v.2.2.29; (Gish and States, 1993; UniProt Consortium, 2019)]. The Uniprot Swiss-Prot Database (downloaded from SwissProt 2018-06-15) was used to obtain protein information and Uniprot Accession codes. Genbank IDs provided by NCBI were used to match CpG background-mRNA overlaps with the annotated mRNA track. Gene ontology terms were paired to Uniprot Accession codes using the Uniprot Swiss-Prot Database (UniProt Consortium, 2019). All GO-MWU inputs are available in the associated Github repository (Venkataraman, 2020).

Once analysis inputs were created, gene ontology terms for each gene were matched with parental terms using default GO-MWU settings. Parental ontology categories with the exact same gene list were combined. Groups were further combined if they shared at least 75% of the same genes. After clustering was complete, a Mann-Whitney U test identified gene ontology categories that were significantly enriched by corresponding CpG loci in genes using the default 10% FDR. Genes with DML were mapped to gene ontology subsets (GO Slim terms) for biological processes to further categorize gene functions.

Results

Water Chemistry

All oysters were initially subjected to acclimation pCO₂ conditions (pCO₂ = 521 ± 32 ppm, Ω_calcite = 2.82 ± 0.13) for 14 days. Following acclimation the treatments were initiated. Oysters in control pCO₂ conditions (pCO₂ = 492 ± 50 μatm; Ω_calcite = 3.01 ± 0.25) experienced low pCO₂ and higher Ω_calcite than those in elevated pCO₂ conditions (pCO₂ = 2550 ± 211 μatm; Ω_calicite = 0.72 ± 0.06) (Table 1).

TABLE 1

Table 1. Summary of water chemistry during the 14-day acclimation period and 28-day experimental exposure.

MBD-BS-Seq

DNA sequencing yielded 280 million DNA sequence reads (NCBI Sequence Read Archive: BioProject accession number PRJNA513384). Of 276 million trimmed paired-end reads, 136 million (49.4%) were mapped to the C. virginica genome, providing an average of 13.6 million reads per sample. Sequencing efforts provided data for 4,304,257 CpG loci (30.7% of 14,458,703 total CpGs in the C. virginica genome) with at least 5× coverage across all samples combined. As expected, the location of CpGs with 5× coverage in the genome differed from the distribution of all CpG motifs (Contingency test; χ² = 1,306,900, df = 6, P-value < 2.2e-16). Of all loci with 5x coverage, 3,255,049 CpGs (75.6%) were found in genic regions in 33,126 out of 38,929 annotated genes in the genome.

The general methylation landscape was defined using all loci with a minimum 5× coverage in each sample. The majority, 3,181,904 (73.9% of MBD-Enriched loci) loci were methylated, with 481,788 (11.2%) sparsely methylated loci and 640,565 (14.9%) unmethylated loci (Figure 1). Median values for global percent methylation and sample methylation varied across genome features (Figure 1). Based on these parameters and data, we calculated that 22% of all CpGs in the gonads (2.7% of total cytosines) had methylation levels greater than 50%. Loci methylation was characterized in relation to putative promoters, UTR, exons, introns, transposable elements, and intergenic regions (Figure 2). Methylated CpGs were found primarily in genic regions, with 2,521,653 loci (79.2%) in 25,496 genes. We rejected the null hypothesis that CpG methylation status was independent of genomic location, as the proportion of methylated CpG loci was different than expected in putative promoters, UTR, exons, introns, transposable elements, and intergenic regions (Contingency test; χ² = 1,311,600, df = 6, P-value < 2.2e-16; Figure 2). There was a larger proportion of methylated loci found in exons compared to all CpGs in the genome (Figure 2). Methylated loci were also found in introns [with 1,448,786 loci (47.3% of methylated loci) vs.1,013,691 CpGs (31.9%) in exons], although this was not higher than expected based on the distribution of all CpGs. Transposable elements contained 755,222 methylated CpGs (23.7%). Putative promoter regions overlapped with 106,111 loci (3.3%), UTR with 128,585 loci (4.0%), and intergenic regions with 660,197 loci (20.7%). There were 372,047 methylated loci (11.7%) that did not overlap with either exons, introns, transposable elements, or promoter regions.

FIGURE 1

Figure 1. Frequency distribution of methylation ratios for CpG loci in C. virginica gonad tissue DNA subjected to MBD enrichment. A total of 4,304,257 CpGs with at least 5× coverage summed across all ten samples were characterized. Loci were considered methylated if they were at least 50% methylated, sparsely methylated loci were 10–50% methylated, and unmethylated loci were 0–10% methylated.

FIGURE 2

Figure 2. Proportion of CpG loci within genomic features. All CpGs are every dinucleotide in the C. virginica genome. Methylated CpGs refers to a dinucleotide with a methylation level of at least 50%.

A total of 37,063 methylation islands were identified in the C. virginica genome (Venkataraman, 2020). Methylation islands contained between 11 and 24,777 methylated CpGs, with a median of 30 methylated CpGs per island. Lengths of methylation islands ranged from 500 to 1,236,482 base pairs, with a median length of 1,024 bp. The majority of methylation islands (36,017; 97.2%) were less than 100,000 bp in length. There were 30,773 (83.0%) methylation islands that overlapped with genic regions.

Differential Methylation Analysis

A total of 598 CpG loci were differentially methylated between oysters exposed to control or high pCO₂, with 51.8% hypermethylated and 48.2% hypomethylated between treatments (Figure 3; Venkataraman, 2020). When considering a PCA using methylation status of all CpG loci with 5× coverage across all samples, the first two principal components explained 29.8% of sample variation (Figure 4A). The first two principal components in a PCA with only differentially methylated loci (DML) explained 57.1% of the variation among treatments (Figure 4B). These DML were distributed throughout the C. virginica genome (Figure 5). The fifth chromosome had the most DML normalized by number of CpGs in the chromosome, and had the most genes; however, this was not the largest chromosome (Figure 5A).

FIGURE 3

Figure 3. Heatmap of DML in C. virginica reproductive tissue created using a euclidean distance matrix. Samples in control pCO₂ conditions are represented by gray, and samples in elevated pCO₂ conditions are represented by a black bar. Loci with higher percent methylation are represented by darker colors. A logistic regression identified 598 DML, defined as individual CpG dinucleotide with at least a 50% methylation change between treatment and control groups, and a q-value < 0.01 based on correction for false discovery rate with the SLIM method. The density of DML at each percent methylation value is represented in the heatmap legend.

FIGURE 4

Figure 4. Principal Components Analysis of (A) all CpG loci with 5x coverage across samples and (B) DML. Methylation status of individual CpG loci explained 29.2% of variation between samples when considering all CpG loci. Methylation status of DML explained 57.1% of sample variation.

FIGURE 5

Figure 5. Distribution of DML among chromosomes and genes. (A) Number of DML normalized by number of CpG in each chromosome (bars) and number of genes (line) in each C. virginica chromosome. (B) Number of genes with various numbers of DML per gene (1–5). Most genes that contained DML only had 1 DML. (C) Proportion of hypermethylated, hypomethylated DML in genes with various numbers of DML per gene (1–5). Mixed refers to a classification of a gene that has both hypermethylated and hypomethylated DML.

Examination of DML within genes revealed that some genes contained multiple DML (Figures 5B,C). Of the 481 genes with DML, the majority only contained one DML (Figure 5B). There were 48 genes with 2 DML, 16 genes with 3 DML, 6 genes with 4 DML and 1 gene with 5 DML (Figure 5B). When multiple DML were found within a gene, they could be methylated in either the same or opposite directions (Figure 5C).

Within the genome, DML were mostly present in genic regions, with 560 DML in 481 genes (368 DML in exons and 192 in introns). In addition, 42 DML were found in putative promoter regions, 27 in UTR, 57 in transposable elements, and 38 in intergenic regions. There were 21 DML located outside of exons, introns, transposable elements, and putative promoters. Additionally, 537 DML were found in methylation islands. The distribution of DML in C. virginica gonad tissue was higher in exons than expected for MBD-enriched CpG loci with minimum 5× coverage across all samples (Contingency test; χ² = 401.09, df = 6, P-value < 2.2e-16; Figure 6). Of the 598 DML, 310 were hypermethylated and 288 were hypomethylated in the high pCO₂ treatment. The number of hyper- and hypomethylated DML was almost evenly split within each genomic feature, with the exception of putative promoter regions that had 44 hypermethylated DML vs. 23 hypomethylated DML. Within a gene, DML did not appear to be concentrated in one particular region. The distribution of hyper- and hypomethylated DML along a gene do not differ from each other (Figure 7).

FIGURE 6

Figure 6. Proportion CpG loci within putative promoters, untranslated regions (UTR), exons, introns, transposable elements, and intergenic regions for MBD-enriched CpGs and differentially methylated loci (DML). The distribution of DML in C. virginica gonad tissue in response to ocean acidification differed from distribution of MBD-enriched loci with 5× coverage across control and treatment samples (Contingency test; χ² = 401.09, df = 6, P-value < 2.2e-16).

FIGURE 7

Figure 7. Distribution of hyper- and hypomethylated DML along a hypothetical gene. The scaled position of a DML within a gene was calculated by dividing the base pair position of the DML by gene length. Counts of hypermethylated DML are plotted above the x-axis, and hypomethylated DML counts are below the x-axis.

The DML were found in genes responsible for various biological processes. However, no gene ontology categories were significantly represented (Figure 8). The majority of genes with DML were involved in protein ubiquitination processes. These genes were not consistently hyper- or hypomethylated. Certain biomineralization genes did contain DML. The gene coding for calmodulin-regulated spectrin-associated protein contained three hypomethylated and one hypermethylated DML. Genes coding for EF-hand protein with calcium-binding domain, calmodulin-binding transcription activator, and calmodulin-lysine N-methyltransferase contained one or two hypermethylated DML.

FIGURE 8

Figure 8. Biological processes represented by all genes used in enrichment background (% Genes) and those with DML (% Genes with DML). Gene ontology categories with similar functions are represented by the same color. Genes may be involved in multiple biological processes. No gene ontologies were significantly enriched.

Discussion

The present study is a general description of DNA methylation in C. virginica, and is one of the first to examine epigenetic responses to ocean acidification in the gonad tissue of a mollusc species. Five hundred ninety-eight differentially methylated loci (DML) were identified in response to the elevated pCO₂ treatments, most of which were in exons. Not only was DNA methylation of C. virginica altered in response to ocean acidification, but changes in gonad methylation indicates potential for these methylation patterns to be inherited by offspring.

Understanding how environmental stressors influence the epigenome is crucial when considering potential acclimatization mechanisms in marine invertebrates. Our finding that high pCO₂ impacts C. virginica DNA methylation adds to a growing body of work about ocean acidification’s impact on marine invertebrate methylomes. The coral species P. damicornis demonstrated an overall increase in DNA methylation when exposed to low pH conditions (7.3–7.6) for 6 weeks, potentially influencing biomineralization (Putnam et al., 2016). Another coral species, S. pistillata, also demonstrated an increase in genome-wide DNA methylation when exposed to low pH conditions for 2 years. Changes in the methylome also modified gene expression and altered pathways involved in cell cycle regulation (Liew et al., 2018b). The present study on an oyster, however, did not observe the overall genome-wide increase in methylation that was reported for corals. Instead, we found subtle, but significant, increases or decreases in percent methylation at several hundred individual CpGs distributed across the genome. As C. virginica and coral species are adapted to different environments and ecological niches, it is possible that species-specific differences in methylation responses contribute to the observed methylation pattern.

The C. virginica methylation landscape suggests a role for methylation in gene activity. Approximately 22% of CpGs in the C. virginica gonad genome were methylated, which is consistent with previous studies of marine invertebrate genomes Gavery and Roberts, 2013; Olson and Roberts, 2014; Hofmann, 2017; Dimond and Roberts, 2020). Methylated loci were concentrated in introns for C. virginica, followed by exons and transposable elements. This location of methylated CpGs in gene bodies is consistent with what has been reported across similar taxa (Roberts and Gavery, 2012; Eirin-Lopez and Putnam, 2018). The concentration of methylated CpGs in gene bodies corresponds with proposed functionality in influencing gene activity (Roberts and Gavery, 2012; Dixon et al., 2014; Liew et al., 2018b). Our study also found methylation in transposable elements, putative promoters and intragenic regions. In plants, transposable element methylation has been shown to modulate the effect of transposable element insertion in genic regions (Hosaka and Kakutani, 2018). It is possible that methylation of transposable elements in C. virginica could also limit the effect of transposable elements. The characterization of methylation islands in the C. virginica genome demonstrates the viability of this descriptive tool for future work examining methylation in mollusc species.

The presence of DML suggests that exposure to experimental ocean acidification conditions elicits an epigenetic response. Many studies have documented changes to oyster protein synthesis, energy production, metabolism, antioxidant responses, and reproduction in response to ocean acidification (Tomanek et al., 2011; Timmins-Schiffman et al., 2014; Dineshram et al., 2016; Boulais et al., 2017; Omoregie et al., 2019). Examination of methylation associated with these physiological responses could identify mechanisms that contribute to these changes. For example, our study found a hypomethylated DML in the heat shock protein 75 kDA gene, and gene expression responses to ocean acidification in C. virginica have found downregulation in a similar molecular chaperone, heat shock protein 70 kDa (Beniash et al., 2010; Ivanina et al., 2014). Other gene expression studies in bivalves have found changes in oxidative stress proteins such as superoxide dismutase, cytochrome c, peroxiredoxin, and NADH dehydrogenase (Chapman et al., 2011; Clark et al., 2013; Goncalves et al., 2016, 2017). Although we did not find any DML in these genes, combined study of DNA methylation and transcription may reveal how changes in gene expression are regulated in response to environmental stressors.

Although DML were found across various genome features, they were mostly in exons and introns. This is consistent with a recent study of C. virginica gill tissue found differentially methylated regions in response to a salinity gradient were primarily in genic regions (Johnson and Kelly, 2019). Interestingly, DML were not found consistently in one particular region of a gene. Similarly, methylated positions in genic regions were evenly distributed after the coral S. pistillata was exposed to low pH (Liew et al., 2018b). Examination of another coral, P. daedalea, in different temperature and salinity conditions found more frequent methylation at 5′ and 3′ ends of genes (Liew et al., 2018a). We also found several genes with multiple DML. These DML were not consistently hyper- or hypomethylated in the same gene. As hyper- and hypomethylation may result in different transcriptional outcomes, future work should examine the role of multiple DML on alternative splicing and gene expression.

The concentration of DML in gene bodies suggests a role for DNA methylation in gene expression and regulation. A majority of genes with DML were involved in protein ubiquitination. Protein ubiquitination is a post-translational protein modification that is involved in protein synthesis and degradation (Peng et al., 2003; Komander, 2009). Previous studies in which oysters were exposed to experimental ocean acidification conditions have demonstrated changes in this pathway. For example, shotgun proteomic characterization of posterior gill lamellae from adult C. gigas exposed to high pCO₂ revealed increased abundance of proteins involved in ubiquitination and decreased protein degradation (Timmins-Schiffman et al., 2014). Elevated pCO₂ levels were also found to upregulate malate dehydrogenase in adult C. virginica mantle tissue (Tomanek et al., 2011). Several genes involved in protein ubiquitination, including those for malate dehydrogenase, ubiquitin-protein ligase, RNA polymerase-associated protein, and DNA damage-binding protein, were significantly hypermethylated in gonad tissue exposed to elevated pCO₂. Hypermethylation of these genes may decrease transcriptional opportunities, thus indicating a critical role in the response to ocean acidification.

Four genes involved in biomineralization contained DML, suggesting these genes can be epigenetically regulated. Upregulation of calcium-binding gene expression has been previously documented in C. virginica (Richards et al., 2018). Since the hypermethylated DML in these genes are typically associated with reduced transcriptional opportunities, it is unclear how methylation changes relate to gene expression for biomineralization genes. Many studies examining ocean acidification-induced carryover effects in bivalves note changes to calcification processes. For example, the Sydney rock oyster (S, glomerata) larvae exhibit faster shell growth in high pCO₂ conditions when parents mature in those same conditions (Parker et al., 2012, 2015). In contrast, larvae from other species found in the North Atlantic such as northern quahog (hard clam; M. mercenaria) and bay scallops (A. irradians) developed slower when parents were reproductively conditioned in low pH conditions (Griffith and Gobler, 2017). There is some evidence to suggest that C. virginica larvae may be more resilient to high pCO₂ conditions than M. mercenaria or A. irradians (Gobler and Talmage, 2014). Differential methylation of biomineralization genes in C. virginica reproductive tissue could be a mechanism to explain when parental experience impacts larval calcification if in fact these DML are inherited.

Although our work documents significant changes to DNA methylation in reproductive tissue after high pCO₂ exposure, this finding may be confounded by secondary effects of gonad maturation. Specimens collected were from mixed populations, and sampled tissue contained both mature and immature gametes. Reproductive tissue likely contained both gametic and somatic cell types. Sex-specific effects have also been documented in response to ocean acidification in mollusc species (Parker et al., 2018; Venkataraman et al., 2019). Lack of a reproductive phenotype precludes any interpretation of how maturation stage or sex can influence changes DNA methylation, as previous work in C. gigas demonstrates these factors as significant influences on baseline methylation patterns (Zhang et al., 2018). Nevertheless, differential methylation in stress response and biomineralization genes suggests that our study does record epigenetic responses to ocean acidification. Future work should pair methylation data with reproductive phenotypes to provide additional information on sex- or stage-specific epigenetic responses to ocean acidification.

Conclusion

Our study found that C. virginica demonstrates a significant epigenetic response to elevated pCO₂ exposure, with 598 DML identified. The concentration of these DML in gene bodies suggests that methylation may be important for transcriptional control in response to environmental stressors. As ocean acidification induced differential methylation in C. virginica gonad tissue, there is a potential for intergenerational epigenetic inheritance, which could control the gene activity of processes such as biomineralization. As carryover effects can persist even when stressors are long-removed (Venkataraman et al., 2019), understanding the mechanisms involved in intergenerational acclimatization is crucial. Future work should focus on methylation patterns in adult C. virginica fully-formed gametes and larvae exposed to various pCO₂ conditions to determine to what degree a difference in methylation influences gene activity and how this might influence phenotypic plasticity.

Data Availability Statement

Raw sequence data is available at the NCBI Sequence Read Archive under BioProject accession number PRJNA513384, with associated metadata and information also available at Woods Hole Open Access Server: https://hdl.handle.net/1912/25138. Associated information for all analyses and supplemental material can be found in the Github repository which is available in an archival format (Venkataraman, 2020; https://doi.org/10.6084/m9.figshare.11923479).

Author Contributions

AD-W, JR, IW, and KL conceived and ran the experiment. YV, SW, and SR performed DNA extractions and methylation analysis. AD-W and KL contributed to analysis. YV and AD-W wrote the initial manuscript draft. All authors reviewed and approved the final manuscript.

Funding

This project was funded by National Science Foundation Biological Oceanography award 1635423 to KL, JR, and BR, and a Hall Conservation Genetics Research Award to YV.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Mackenzie Gavery and our two reviewers for their insightful feedback on the manuscript. This work was facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington.

Footnotes

1. ^ https://github.com/z0on/GO_MWU

References

Akalin, A., Kormaksson, M., Li, S., Garrett-Bakelman, F. E., Figueroa, M. E., Melnick, A., et al. (2012). methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 13:R87.

Google Scholar

Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data [Online]. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

Google Scholar

Bao, W., Kojima, K. K., and Kohany, O. (2015). Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob. DNA 6:11.

Google Scholar

Barton, A., Hales, B., Waldbusser, G. G., Langdon, C., and Feely, R. A. (2012). The Pacific oyster, Crassostrea gigas, shows negative correlation to naturally elevated carbon dioxide levels: implications for near-term ocean acidification effects. Limnol. Oceanogr. 57, 698–710. doi: 10.4319/lo.2012.57.3.0698

CrossRef Full Text | Google Scholar

Beniash, E., Ivanina, A., Lieb, N. S., Kurochkin, I., and Sokolova, I. M. (2010). Elevated level of carbon dioxide affects metabolism and shell formation in oysters Crassostrea virginica. Mar. Ecol. Prog. Ser. 419, 95–108. doi: 10.3354/meps08841

CrossRef Full Text | Google Scholar

Bird, A. (2002). DNA methylation patterns and epigenetic memory. Genes Dev. 16, 6–21. doi: 10.1101/gad.947102

CrossRef Full Text | Google Scholar

Bossdorf, O., Richards, C. L., and Pigliucci, M. (2008). Epigenetics for ecologists. Ecol. Lett. 11, 106–115.

Google Scholar

Boulais, M., Chenevert, K. J., Demey, A. T., Darrow, E. S., Robison, M. R., Roberts, J. P., et al. (2017). Oyster reproduction is compromised by acidification experienced seasonally in coastal regions. Sci. Rep. 7:13276.

Google Scholar

Byrne, M., Foo, S. A., Ross, P. M., and Putnam, H. M. (2019). Limitations of cross and multigenerational plasticity for marine invertebrates faced with global climate change. Glob. Chang. Biol. 26, 80–102. doi: 10.1111/gcb.14882

CrossRef Full Text | Google Scholar

Chapman, R. W., Mancia, A., Beal, M., Veloso, A., Rathburn, C., Blair, A., et al. (2011). The transcriptomic responses of the eastern oyster, Crassostrea virginica, to environmental conditions. Mol. Ecol. 20, 1431–1449. doi: 10.1111/j.1365-294x.2011.05018.x

CrossRef Full Text | Google Scholar

Clark, M. S., Thorne, M. A. S., Amaral, A., Vieira, F., Batista, F. M., Reis, J., et al. (2013). Identification of molecular and physiological responses to chronic environmental challenge in an invasive species: the Pacific oyster, Crassostrea gigas. Ecol. Evol. 3, 3283–3297.

Google Scholar

Deans, C., and Maggert, K. A. (2015). What do you mean, “Epigenetic”? Genetics 199, 887–896. doi: 10.1534/genetics.114.173492

CrossRef Full Text | Google Scholar

Dickson, A. G. (1990). Standard potential of the reaction : AgCl(s) + 12H2(g) = Ag(s) + HCl(aq), and the standard acidity constant of the ion HSO4- in synthetic sea water from 273.15 to 318.15 K. J. Chem. Thermodyn. 22, 113–127. doi: 10.1016/0021-9614(90)90074-z

CrossRef Full Text | Google Scholar

Dimond, J. L., and Roberts, S. B. (2016). Germline DNA methylation in reef corals: patterns and potential roles in response to environmental change. Mol. Ecol. 25, 1895–1904. doi: 10.1111/mec.13414

CrossRef Full Text | Google Scholar

Dimond, J. L., and Roberts, S. B. (2020). Convergence of DNA methylation profiles of the reef coral Porites astreoides in a novel environment. Front. Mar. Sci. 6:792. doi: 10.3389/fmars.2019.00792

CrossRef Full Text | Google Scholar

Dineshram, R., Chandramouli, K., Ko, G. W. K., Zhang, H., Qian, P.-Y., Ravasi, T., et al. (2016). Quantitative analysis of oyster larval proteome provides new insights into the effects of multiple climate change stressors. Glob. Chang. Biol. 22, 2054–2068. doi: 10.1111/gcb.13249

CrossRef Full Text | Google Scholar

Dixon, G. B., Bay, L. K., and Matz, M. V. (2014). Bimodal signatures of germline methylation are linked with gene expression plasticity in the coral Acropora millepora. BMC Genomics 15:1109. doi: 10.1186/1471-2164-15-1109

CrossRef Full Text | Google Scholar

Eirin-Lopez, J. M., and Putnam, H. M. (2018). Marine environmental epigenetics. Ann. Rev. Mar. Sci. 11, 335–368. doi: 10.1146/annurev-marine-010318-095114

CrossRef Full Text | Google Scholar

Ekstrom, J. A., Suatoni, L., Cooley, S. R., Pendleton, L. H., Waldbusser, G. G., Cinner, J. E., et al. (2015). Vulnerability and adaptation of US shellfisheries to ocean acidification. Nat. Clim. Chang. 5, 207–214. doi: 10.1038/nclimate2508

CrossRef Full Text | Google Scholar

Feely, R. A., Alin, S. R., Newton, J., Sabine, C. L., Warner, M., Devol, A., et al. (2010). The combined effects of ocean acidification, mixing, and respiration on pH and carbonate saturation in an urbanized estuary. Estuar. Coast. Shelf Sci. 88, 442–449. doi: 10.1016/j.ecss.2010.05.004

CrossRef Full Text | Google Scholar

Gatzmann, F., Falckenhayn, C., Gutekunst, J., Hanna, K., Raddatz, G., Carneiro, V. C., et al. (2018). The methylome of the marbled crayfish links gene body methylation to stable expression of poorly accessible genes. Epigenet. Chromatin 11:57.

Google Scholar

Gavery, M. R., and Roberts, S. B. (2013). Predominant intragenic methylation is associated with gene expression characteristics in a bivalve mollusc. PeerJ 1:e215. doi: 10.7717/peerj.215

CrossRef Full Text | Google Scholar

Gazeau, F., Quiblier, C., Jansen, J. M., Gattuso, J.-P., Middelburg, J. J., and Heip, C. H. R. (2007). Impact of elevated CO2 on shellfish calcification. Geophys. Res. Lett. 34:L07603.

Google Scholar

Gish, W., and States, D. J. (1993). Identification of protein coding regions by database similarity search. Nat. Genet. 3, 266–272. doi: 10.1038/ng0393-266

CrossRef Full Text | Google Scholar

Gobler, C. J., and Talmage, S. C. (2014). Physiological response and resilience of early life-stage Eastern oysters (Crassostrea virginica) to past, present and future ocean acidification. Conserv. Physiol. 2:cou004. doi: 10.1093/conphys/cou004

CrossRef Full Text | Google Scholar

Gómez-Chiarri, M., Warren, W. C., Guo, X., and Proestou, D. (2015). Developing tools for the study of molluscan immunity: the sequencing of the genome of the eastern oyster, Crassostrea virginica. Fish Shellfish Immunol. 46, 2–4. doi: 10.1016/j.fsi.2015.05.004

CrossRef Full Text | Google Scholar

Goncalves, P., Anderson, K., Thompson, E. L., Melwani, A., Parker, L. M., Ross, P. M., et al. (2016). Rapid transcriptional acclimation following transgenerational exposure of oysters to ocean acidification. Mol. Ecol. 25, 4836–4849. doi: 10.1111/mec.13808

CrossRef Full Text | Google Scholar

Goncalves, P., Thompson, E. L., and Raftos, D. A. (2017). Contrasting impacts of ocean acidification and warming on the molecular responses of CO2-resilient oysters. BMC Genomics 18:431. doi: 10.1186/s12864-017-3818-z

CrossRef Full Text | Google Scholar

Griffith, A. W., and Gobler, C. J. (2017). Transgenerational exposure of North Atlantic bivalves to ocean acidification renders offspring more vulnerable to low pH and additional stressors. Sci. Rep. 7:11394.

Google Scholar

Helm, M. M., and Bourne, N. (2004). Hatchery Culture of Bivalves: A Practical Manual. Rome: Food & Agriculture Org.

Google Scholar

Hofmann, G. E. (2017). Ecological epigenetics in marine metazoans. Front. Mar. Sci. 4:4. doi: 10.3389/fmars.2017.00004

CrossRef Full Text | Google Scholar

Hosaka, A., and Kakutani, T. (2018). Transposable elements, genome evolution and transgenerational epigenetic variation. Curr. Opin. Genet. Dev. 49, 43–48. doi: 10.1016/j.gde.2018.02.012

CrossRef Full Text | Google Scholar

IPCC (2019). “Summary for Policymakers,” in IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, eds H.-O. Pörtner, D. C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, et al. (Geneva: IPCC).

Google Scholar

Ivanina, A. V., Hawkins, C., and Sokolova, I. M. (2014). Immunomodulation by the interactive effects of cadmium and hypercapnia in marine bivalves Crassostrea virginica and Mercenaria mercenaria. Fish Shellfish Immunol. 37, 299–312. doi: 10.1016/j.fsi.2014.02.016

CrossRef Full Text | Google Scholar

Jeong, H., Wu, X., Smith, B., and Yi, S. V. (2018). Genomic Landscape of Methylation Islands in Hymenopteran Insects. Genome Biol. Evol. 10, 2766–2776. doi: 10.1093/gbe/evy203

CrossRef Full Text | Google Scholar

Johnson, K. M., and Kelly, M. W. (2019). Population epigenetic divergence exceeds genetic divergence in the Eastern oyster Crassostrea virginica in the Northern Gulf of Mexico. Evol. Appl. 00, 1–15. doi: 10.1111/eva.12912

CrossRef Full Text | Google Scholar

Komander, D. (2009). The emerging complexity of protein ubiquitination. Biochem. Soc. Trans. 37, 937–953. doi: 10.1042/bst0370937

CrossRef Full Text | Google Scholar

Krueger, F., and Andrews, S. R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572. doi: 10.1093/bioinformatics/btr167

CrossRef Full Text | Google Scholar

Kurihara, H., Kato, S., and Ishimatsu, A. (2007). Effects of increased seawater pCO2 on early development of the oyster Crassostrea gigas. Aquat. Biol. 1, 91–98. doi: 10.3354/ab00009

CrossRef Full Text | Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923

CrossRef Full Text | Google Scholar

Lee, K., Tae-Wook, K., Byrne, R., Millero, F. J., Feely, R. A., and Liu, Y.-M. (2010). The universal ratio of boron to chlorinity for the North Pacific and North Atlantic oceans. Geochim. Cosmochim. Acta 74, 1801–1811. doi: 10.1016/j.gca.2009.12.027

CrossRef Full Text | Google Scholar

Lewis, E., and Wallace, D. W. (1998). R: Program developed for CO2 system calculations ORNL/CDIAC-105. Oak Ridge: US Department of Energy.

Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

CrossRef Full Text | Google Scholar

Liew, Y. J., Howells, E. J., Wang, X., Michell, C. T., Burt, J. A., Idaghdour, Y., et al. (2018a). Intergenerational epigenetic inheritance in reef-building corals. bioRxiv [preprint]. doi: 10.1101/269076

CrossRef Full Text | Google Scholar

Liew, Y. J., Zoccola, D., Li, Y., Tambutté, E., Venn, A. A., and Michell, C. T. (2018b). Epigenome-associated phenotypic acclimatization to ocean acidification in a reef-building coral. Sci. Adv. 4:eaar8028. doi: 10.1126/sciadv.aar8028

CrossRef Full Text | Google Scholar

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12.

Google Scholar

Olson, C. E., and Roberts, S. B. (2014). Genome-wide profiling of DNA methylation and gene expression in Crassostrea gigas male gametes. Front. Physiol. 5:224. doi: 10.3389/fphys.2014.00224

CrossRef Full Text | Google Scholar

Omoregie, E., Mwatilifange, N. S. I., and Liswaniso, G. (2019). Futuristic ocean acidification levels reduce growth and reproductive viability in the pacific oyster (Crassostrea gigas). J. Appl. Sci. Environ. Manage. 23, 1747–1754.

Google Scholar

Parker, L. M., O’Connor, W. A., Byrne, M., Coleman, R. A., Virtue, P., Dove, M., et al. (2017). Adult exposure to ocean acidification is maladaptive for larvae of the Sydney rock oyster Saccostrea glomerata in the presence of multiple stressors. Biol. Lett. 13:20160798. doi: 10.1098/rsbl.2016.0798

CrossRef Full Text | Google Scholar

Parker, L. M., O’Connor, W. A., Byrne, M., Dove, M., Coleman, R. A., Pörtner, H.-O., et al. (2018). Ocean acidification but not warming alters sex determination in the Sydney rock oyster, Saccostrea glomerata. Proc. R. Soc. B 285:20172869. doi: 10.1098/rspb.2017.2869

CrossRef Full Text | Google Scholar

Parker, L. M., O’Connor, W. A., Raftos, D. A., Pörtner, H.-O., and Ross, P. M. (2015). Persistence of positive carryover effects in the Oyster, Saccostrea glomerata, following transgenerational exposure to ocean acidification. PLoS One 10:e0132276. doi: 10.1371/journal.pone.0132276

CrossRef Full Text | Google Scholar

Parker, L. M., Ross, P. M., O’Connor, W. A., Borysko, L., Raftos, D. A., and Pörtner, H.-O. (2012). Adult exposure influences offspring response to ocean acidification in oysters. Glob. Chang. Biol. 18, 82–92. doi: 10.1111/j.1365-2486.2011.02520.x

CrossRef Full Text | Google Scholar

Parker, L. M., Ross, P. M., O’Connor, W. A., Pörtner, H. O., Scanes, E., and Wright, J. M. (2013). Predicting the response of molluscs to the impact of ocean acidification. Biology 2, 651–692. doi: 10.3390/biology2020651

CrossRef Full Text | Google Scholar

Peng, J., Schwartz, D., Elias, J. E., Thoreen, C. C., Cheng, D., Marsischky, G., et al. (2003). A proteomics approach to understanding protein ubiquitination. Nat. Biotechnol. 21, 921–926. doi: 10.1038/nbt849

CrossRef Full Text | Google Scholar

Putnam, H. M., Davidson, J. M., and Gates, R. D. (2016). Ocean acidification influences host DNA methylation and phenotypic plasticity in environmentally susceptible corals. Evol. Appl. 9, 1165–1178. doi: 10.1111/eva.12408

CrossRef Full Text | Google Scholar

Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. doi: 10.1093/bioinformatics/btq033

CrossRef Full Text | Google Scholar

Richards, M., Xu, W., Mallozzi, A., Errera, R. M., and Supan, J. (2018). Production of calcium-binding proteins in Crassostrea virginica in response to increased environmental CO2 concentration. Front. Mar. Sci. 5:203. doi: 10.3389/fmars.2018.00203

CrossRef Full Text | Google Scholar

Ries, J. B. (2011). A physicochemical framework for interpreting the biological calcification response to CO2-induced ocean acidification. Geochim. Cosmochim. Acta 75, 4053–4064. doi: 10.1016/j.gca.2011.04.025

CrossRef Full Text | Google Scholar

Roberts, S. B., and Gavery, M. R. (2012). Is There a Relationship between DNA Methylation and Phenotypic Plasticity in Invertebrates? Front. Physiol. 2:116. doi: 10.3389/fphys.2011.00116

CrossRef Full Text | Google Scholar

Rondon, R., Grunau, C., Fallet, M., Charlemagne, N., Sussarellu, R., Chaparro, C., et al. (2017). Effects of a parental exposure to diuron on Pacific oyster spat methylome. Environ. Epigenet. 3:dvx004. doi: 10.1093/eep/dvx004

CrossRef Full Text | Google Scholar

Ross, P. M., Parker, L., and Byrne, M. (2016). Transgenerational responses of molluscs and echinoderms to changing ocean conditions. ICES J. Mar. Sci. 73, 537–549. doi: 10.1093/icesjms/fsv254

CrossRef Full Text | Google Scholar

Roy, R. N., Roy, L. N., Vogel, K. M., Porter-Moore, C., Pearson, T., Good, C. E., et al. (1993). The dissociation constants of carbonic acid in seawater at salinities 5 to 45 and temperatures 0 to 45°C. Mar. Chem. 44, 249–267. doi: 10.1016/0304-4203(93)90207-5

CrossRef Full Text | Google Scholar

Smit, A. F. A., Hubley, R., and Green, P. (2013). RepeatMasker Open-4.0. Available online at: http://www.repeatmasker.org

Google Scholar

Strader, M. E., Wong, J. M., Kozal, L. C., Leach, T. S., and Hofmann, G. E. (2019). Parental environments alter DNA methylation in offspring of the purple sea urchin, Strongylocentrotus purpuratus. J. Exp. Mar. Biol. Ecol. 517, 54–64. doi: 10.1016/j.jembe.2019.03.002

CrossRef Full Text | Google Scholar

Suzuki, M. M., and Bird, A. (2008). DNA methylation landscapes: provocative insights from epigenomics. Nat. Rev. Genet. 9, 465–476. doi: 10.1038/nrg2341

CrossRef Full Text | Google Scholar

Timmins-Schiffman, E., Coffey, W. D., Hua, W., Nunn, B. L., Dickinson, G. H., and Roberts, S. B. (2014). Shotgun proteomics reveals physiological response to ocean acidification in Crassostrea gigas. BMC Genomics 15:951. doi: 10.1186/1471-2164-15-951

CrossRef Full Text | Google Scholar

Tomanek, L., Zuzow, M. J., Ivanina, A. V., Beniash, E., and Sokolova, I. M. (2011). Proteomic response to elevated PCO2 level in eastern oysters, Crassostrea virginica: evidence for oxidative stress. J. Exp. Biol. 214, 1836–1844. doi: 10.1242/jeb.055475

CrossRef Full Text | Google Scholar

UniProt Consortium (2019). UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 47, D506–D515.

Google Scholar

Van Heuven, S., Pierrot, D., Rae, J. W. B., Lewis, E., and Wallace, D. W. R. (2011). MATLAB program developed for CO2 system calculations. ORNL/CDIAC-105b. Oak Ridge: US Department of Energy.

Google Scholar

Venkataraman, Y. R. (2020). Eastern oyster (Crassostrea virginica) gonad DNA methylation data and analysis. Dataset. doi: 10.6084/m9.figshare.11923479.v3

CrossRef Full Text | Google Scholar

Venkataraman, Y. R., Spencer, L. H., and Roberts, S. B. (2019). Larval response to parental low ph exposure in the Pacific Oyster Crassostrea gigas. J. Shellfish Res. 38:743. doi: 10.2983/035.038.0325

CrossRef Full Text | Google Scholar

Waldbusser, G. G., Hales, B., Langdon, C. J., Haley, B. A., Schrader, P., Brunner, E. L., et al. (2014). Saturation-state sensitivity of marine bivalve larvae to ocean acidification. Nat. Clim. Chang. 5:273. doi: 10.1038/nclimate2479

CrossRef Full Text | Google Scholar

Wang, H.-Q., Tuominen, L. K., and Tsai, C.-J. (2011). SLIM: a sliding linear model for estimating the proportion of true null hypotheses in datasets with dependence structures. Bioinformatics 27, 225–231. doi: 10.1093/bioinformatics/btq650

CrossRef Full Text | Google Scholar

Whitlock, M., and Schluter, D. (2014). The Analysis of Biological Data. Chennai: Roberts & Company.

Google Scholar

Wright, R. M., Aglyamova, G. V., Meyer, E., and Matz, M. V. (2015). Gene expression associated with white syndromes in a reef building coral, Acropora hyacinthus. BMC Genomics 16:371. doi: 10.1186/s12864-015-1540-2

CrossRef Full Text | Google Scholar

Zhang, X., Li, Q., Kong, L., and Yu, H. (2018). DNA methylation frequency and epigenetic variability of the Pacific oyster Crassostrea gigas in relation to the gametogenesis. Fish. Sci. 84, 789–797. doi: 10.1007/s12562-018-1214-5

CrossRef Full Text | Google Scholar

Zhao, L., Liu, L., Liu, B., Liang, J., Lu, Y., and Yang, F. (2019). Antioxidant responses to seawater acidification in an invasive fouling mussel are alleviated by transgenerational acclimation. Aquat. Toxicol. 217:105331. doi: 10.1016/j.aquatox.2019.105331

CrossRef Full Text | Google Scholar

Zhao, L., Yang, F., Milano, S., Han, T., Walliser, E. O., and Schöne, B. R. (2018). Transgenerational acclimation to seawater acidification in the Manila clam Ruditapes philippinarum: Preferential uptake of metabolic carbon. Sci. Total Environ. 627, 95–103. doi: 10.1016/j.scitotenv.2018.01.225

CrossRef Full Text | Google Scholar