Ocean Acidification Induces Subtle Shifts in Gene Expression and DNA Methylation in Mantle Tissue of the Eastern Oyster (Crassostrea virginica)

Early evidence suggests that DNA methylation can mediate phenotypic responses of marine calcifying species to ocean acidification (OA). Few studies, however, have explicitly studied DNA methylation in calcifying tissues through time. Here, we examined the phenotypic and molecular responses in the extrapallial fluid and mantle (fluid and tissue at the calcification site) in adult eastern oyster (Crassostrea virginica) exposed to experimental OA over 80 days. Oysters were reared under three experimental pCO2 treatments (“control,” 580 μatm; “moderate OA,” 1,000 μatm; “high OA,” 2,800 μatm) and sampled at 6 time points (24 h−80 days). We found that high OA initially induced an increase in the pH of the extrapallial fluid (pHEPF) relative to the external seawater that peaked at day 9, but then diminished over time. Calcification rates were significantly lower in the high OA treatment compared to the other treatments. To explore how oysters regulate their extrapallial fluid, gene expression and DNA methylation were examined in the mantle-edge tissue of oysters from days 9 and 80 in the control and high OA treatments. Mantle tissue mounted a significant global molecular response (both in the transcriptome and methylome) to OA that shifted through time. Although we did not find individual genes that were significantly differentially expressed under OA, the pHEPF was significantly correlated with the eigengene expression of several co-expressed gene clusters. A small number of OA-induced differentially methylated loci were discovered, which corresponded with a weak association between OA-induced changes in genome-wide gene body DNA methylation and gene expression. Gene body methylation, however, was not significantly correlated with the eigengene expression of pHEPF-correlated gene clusters. These results suggest that OA induces a subtle response in a large number of genes in C. virginica, but also indicate that plasticity at the molecular level may be limited. Our study highlights the need to reassess our understanding of tissue-specific molecular responses in marine calcifiers, as well as the role of DNA methylation and gene expression in mediating physiological and biomineralization responses to OA.

Early evidence suggests that DNA methylation can mediate phenotypic responses of marine calcifying species to ocean acidification (OA). Few studies, however, have explicitly studied DNA methylation in calcifying tissues through time. Here, we examined the phenotypic and molecular responses in the extrapallial fluid and mantle (fluid and tissue at the calcification site) in adult eastern oyster (Crassostrea virginica) exposed to experimental OA over 80 days. Oysters were reared under three experimental pCO 2 treatments ("control," 580 µatm; "moderate OA," 1,000 µatm; "high OA," 2,800 µatm) and sampled at 6 time points (24 h−80 days). We found that high OA initially induced an increase in the pH of the extrapallial fluid (pH EPF ) relative to the external seawater that peaked at day 9, but then diminished over time. Calcification rates were significantly lower in the high OA treatment compared to the other treatments. To explore how oysters regulate their extrapallial fluid, gene expression and DNA methylation were examined in the mantle-edge tissue of oysters from days 9 and 80 in the control and high OA treatments. Mantle tissue mounted a significant global molecular response (both in the transcriptome and methylome) to OA that shifted through time. Although we did not find individual genes that were significantly differentially expressed under OA, the pH EPF was significantly correlated with the eigengene expression of several co-expressed gene clusters. A small number of OA-induced differentially methylated loci were discovered, which corresponded with a weak association between OA-induced changes in genome-wide gene body DNA methylation and gene expression. Gene body methylation, however, was not significantly correlated with the eigengene expression of pH EPF -correlated gene clusters. These results suggest that OA induces a subtle response in a large number of genes in C. virginica, but also indicate that plasticity at the molecular level may be limited. Our study highlights the need to reassess our understanding of tissue-specific molecular responses in marine calcifiers, as well as the role of DNA methylation and gene expression in mediating physiological and biomineralization responses to OA.
Keywords: DNA methylation, gene expression, ocean acidification, acclimatization, Crassostrea virginica, oyster, extrapallial fluid INTRODUCTION Ocean acidification (OA), a decrease in seawater pH due to the uptake of anthropogenic CO 2 , is expected to have substantial effects on marine species and ecosystems in the near future (Orr et al., 2005;Guinotte and Fabry, 2008;Doney et al., 2009). Specifically, OA is driving a shift in the carbonate system equilibrium, resulting in decreased availability of carbonate ions and lower calcium carbonate saturation state (Feely et al., 2004;Orr et al., 2005). This may be particularly problematic for marine calcifying species, which build their shells and skeletons from calcium and carbonate ions. Substantial effort has been invested into evaluating the shortand long-term consequences of OA across taxa and life history stages. Prolonged exposure to OA tends to have negative effects on calcification, metabolism, and growth, as observed in corals (e.g., Anthony et al., 2008), pteropods (e.g., Bednaršek et al., 2012), gastropods (e.g., Melatunan et al., 2013), and bivalves (e.g., Talmage and Gobler, 2010). However, these negative effects are not universal, varying in direction (negative or positive) and intensity depending on taxon, severity of OA, and duration of exposure (Ries et al., 2009;Kroeker et al., 2010). Observed nonnegative effects may be the result of individual resilience or the capacity of certain species to acclimatize to OA.
Under ambient conditions, adult bivalves tend to maintain their extra-pallial fluid pH (pH EPF ) below that of seawater pH (pH seawater ) (Crenshaw, 1972;Sutton et al., 2018), which has been attributed to the buildup of dissolved CO 2 when the organism's valves are closed (Crenshaw, 1972;Michaelidis et al., 2005;Sutton et al., 2018;Liu et al., 2020). While under OA conditions, some bivalves will actually increase their pH EPF above that of pH seawater (e.g., Liu et al., 2020), potentially enabling them to continue calcifying even under highly acidic scenarios. Only recently have investigators studied the molecular mechanisms that regulate calcification within the mantle tissue (see reviews by Clark, 2020;Rajan and Vengatesen, 2020).
Transcriptomic studies have proven to be a powerful approach for identifying genes and understanding pathways that shape organismal the response to OA, and have been used to help elucidate the potential trade-offs that result from acclimatization (Evans and Hofmann, 2012), including genes associated with biomineralization, acid-base regulation, and metabolic function (Evans et al., 2013;Davies et al., 2016;Li et al., 2016;Goncalves et al., 2017;Wong et al., 2018;Griffiths et al., 2019;Strader et al., 2020). From this literature it has been hypothesized that there is a trade-off between maintaining calcification vs. other core functions (Wood et al., 2008), whereby increased or sustained calcification under continued exposure may become too costly over time. At the transcriptomic level, this cost can manifest as shifts in the expression of both calcification and metabolic genes over time (Li et al., 2016), suggesting that certain acclimatization responses may be unsustainable during extended exposure to OA. The precise mechanisms that mediate these changes in gene expression remain underexplored, but recent studies indicate that DNA methylation, an epigenetic modification, may be an important regulator of these responses (Liew et al., 2018;Bogan et al., 2020).
DNA methylation is an important mediator of gene regulation and physiological response across a diverse range of taxa (see reviews Zilberman et al., 2007;Zemach et al., 2010;Schübeler, 2015;Zhang et al., 2018). DNA methylation refers to the attachment of a methyl group to cytosine, and is most frequently studied in the context of the cytosine-guanine (CpG) motif. These modifications occur on top of the genome without altering the underlying DNA sequence (Richards, 2006;Bonduriansky and Day, 2018) and may be stable enough to be passed from parent to offspring (Jablonka and Lamb, 2002). Previous work has shown they can have an important role in determining an organism's phenotype across a range of species, including fur color and obesity in mice (Dolinoy et al., 2006), rat pup behavior as determined by parental care (Anier et al., 2014), response to phosphate starvation in Arabidopsis (Yong-Villalobos et al., 2015), and caste phenotypes in social insects (Lyko et al., 2010). Importantly, the precise manner in which DNA methylation serves to regulate phenotype is taxon-specific.
In invertebrates, association studies have begun to elucidate the relationship between DNA methylation, gene expression, and phenotype (e.g., Gavery and Roberts, 2013;Dixon et al., 2018;Liew et al., 2018). In basal invertebrates (e.g., corals and oysters), genomes are typically only sparsely methylated compared to vertebrates, and DNA methylation tends to be concentrated in CpGs within gene bodies (i.e., exons and introns). In this group, gene body methylation was positively correlated with gene expression and negatively correlated with the variation in expression among individuals (Gavery and Roberts, 2013;Liew et al., 2018). Based on these patterns, DNA methylation may be acting to regulate not only the level of expression, but also the amount of transcriptional noise among individuals (Huh et al., 2013). Furthermore, recent evidence suggests that unlike vertebrates, non-deuterostome invertebrate DNA methylation is not necessarily erased during embryogenesis, opening up the possibility for transmission of marks between generations (Xu et al., 2019).
Understanding the role of environment-induced DNA methylation in regulating gene expression and phenotypic plasticity is of particular interest in the context of marine systems and global change (see reviews Hofmann, 2017;Eirin-Lopez and Putnam, 2019). In corals, early evidence indicates that DNA methylation can be both sensitive to environmental change, including OA in some species (Putnam et al., 2016;Liew et al., 2018), and associated with gene expression (Dixon et al., 2018) suggesting a potential regulatory role of DNA methylation in coral response to OA. Liew et al. (2018) found OA-sensitive DNA methylation also corresponded with biomineralization-related traits in corals, supporting the idea that DNA methylation may have an important role in phenotypic plasticity and acclimatization to OA. Although these studies examined OA reponses in many different tissues, none of these studies quantified methylation responses to OA at base-pair resolution in calcifying tissue.
Here, we examine the molecular basis of eastern oyster (Crassostrea virginica) pH EPF regulation and calcification response to OA. Crassostrea virginica is a coastal marine bivalve that provides critical ecosystem services in the form of habitat structure and water filtration, and also serves as an economically important food source (Ekstrom et al., 2015;Gómez-Chiarri et al., 2015). Its genome has been recently sequenced (Gómez-Chiarri et al., 2015) and the species has been the subject of extensive OA studies. Consistent with other bivalve species, C. virginica at all life stages exhibit a largely negative calcification response to OA (e.g., Miller et al., 2009;Ries et al., 2009;Beniash et al., 2010;Waldbusser et al., 2011b;Dickinson et al., 2012), while metabolic and physiological responses tend to be more varied and dependent on life stage, experimental condition, and duration of exposure (Beniash et al., 2010;Dickinson et al., 2012;Ivanina et al., 2013;Matoo et al., 2013;Gobler and Talmage, 2014). Like other bivalve species, the chemical composition of C. virginica EPF is thought to be regulated via the active exchange of ions and other constituents through the mantle epithelium to promote calcification (Crenshaw and Neff, 1969;Crenshaw, 1972). A recent study by Liu et al. (2020) found that C. virginica had some ability to moderate their pH EPF in response to OA, however, the mechanisms that regulate this response have not been explored.
The primary goal of the present study was to investigate the molecular mechanisms underlying the biomineralization response of C. virginica to short-and long-term OA exposure. Specifically, we studied the role of DNA methylation in physiological and gene expression responses in C. virginica mantle-edge tissue. We performed a controlled laboratory experiment in which oysters were exposed to three levels of CO 2 -induced OA (control, moderate, and high) for 80 days and monitored the oysters' calcification rate and pH EPF . Based on patterns we observed in the pH EPF response, we examined gene expression and DNA methylation responses in mantleedge tissue at two time points (days 9 and 80) for oysters exposed to control and high OA treatments. We predicted that the oysters mitigate the effects of OA by regulating pH EPF relative to pH seawater and that this process is associated with a molecular response in the oyster transcriptome and methylome. Specifically, we hypothesized that in C. virginica (i) OA will drive differential responses in both the transcriptome and methylome and up-regulation of pH EPF relative to the environment, (ii) pH EPF response will be coupled with the up-regulation of biomineralization related genes, (iii) gene expression will be associated with gene body DNA methylation, and (iv) differential expression will correspond to significant shifts in gene body methylation. By establishing the first integrated dataset on EPF carbonate chemistry, gene expression, and DNA methylation in a mollusk exposed to variable degrees of OA, our study provides new insights into the responses of bivalves to OA over short and long timescales.

Oyster Collection and Preparation
Adult C. virginica were collected from three intertidal sites within Plum Island Sound, Massachusetts, USA (Site 1, 42.751636, −70.837023;Site 2, 42.725186, −70.855022;Site 3, 42.681764, −70.813498) in late April 2017. These sites are all within 8 km of each other and the oysters from these sites are not distinct genetic populations. However, to compensate for any site effects, site was included as a random effect in applicable statistical models (see statistical analysis section for details).
In the lab, oysters were cleaned of epibionts and shell ports were installed over 8 days while being maintained in 50-L flowthrough tanks. To measure the carbonate chemistry of the EPF, a 2 mm hole was drilled into the right valve ∼2 cm from the hinge to expose the EPF cavity, without damaging the underlying mantle tissue. The drilled hole was gently rinsed with filtered seawater and patted dry. The barbed end of a nylon luer-lock coupling (McMaster-Carr 51525K123) was trimmed to the length of the shell thickness, inserted into the hole and sealed in place using marine-safe cyanoacrylate (Starbond EM-2000 CA USA). The coupling was then plugged with a nylon cap socket (McMaster-Carr 51525K315) to prevent exchange between seawater and the EPF via the port and limit the potential for increased dissolution near the site of the implant. Oysters were allowed to recover for 4 days, after which they were transferred to the experimental tank array (described below) for acclimation and exposure. During the transfer, oysters were assigned one of three pCO 2 treatments: "control" (ca. 580 µatm)-corresponding to near-present-day conditions in the estuary where they were collected; "moderate OA" (ca. 1,000 µatm)-corresponding to end-of-century predictions (Intergovernmental Panel on Climate Change, 2007); and "high OA" (ca. 2,800 µatm)-corresponding to seawater conditions that are undersaturated with respect to calcite ( calcite < 1). Importantly, high OA events are already regularly experienced in the coastal marine systems from which these oysters were collected, owing to daily tidal cycles, input of meteoric waters, and seasonal reoxidation of organic matter ( Supplementary Figure 1.1). This trend is underappreciated in many coastal systems, where local processes are increasing the acidity of the seawater relative to the global average (Reum et al., 2014). The high OA treatment was also used to increase the inferential and statistical power of the study design (Whitlock and Schluter, 2014).

Experimental Design and Water Chemistry
The experiment was conducted at Northeastern University's Marine Science Center, using a flow-through seawater system that draws water from Broad Sound in Nahant, Massachusetts (42.416884, −70.907564). Each pCO 2 treatment was replicated over two blocks, each containing three 42-L acrylic aquaria (total of 6 aquaria per treatment). Tanks within a block were connected via the recirculation and filtration system, while gas mixture, temperature, and incoming filtered seawater were regulated independently for each tank. A blocked randomized design was used to ensure an equal number of oysters from each location were distributed across all treatments and that there were no significant differences in initial oyster size across treatment (mean ± SD shell length: 9.70 ± 2.4 cm). Twelve oysters were placed into each aquarium, four from each location, resulting in 24 oysters per treatment per location (total of 72 oysters per treatment). Oysters were acclimated for 33 days in control pCO 2 conditions at a temperature of 17 • C. This extended acclimation period was used to (i) provide time for any acute field-induced changes due to local acidification events to reset, and (ii) ensure that the oysters had sufficient time to recover from the stress of being transported to the laboratory. Following acclimation, the moderate and high OA treatments were ramped up to target pCO 2 over a 12-h period. This 12-h period was the minimum amount of time needed for the experimental OA system to equilibrate with the new pCO 2 conditions and was chosen to simulate the rapid change in OA observed from the field (Supplementary Figure 1.1). Treatment conditions were maintained for 80 days.
Target pCO 2 gases were formulated by mixing compressed CO 2 with compressed CO 2 -free air (control treatment) or with compressed air (moderate and high OA treatments) using solenoid-valve-controlled mass flow controllers at flow-rates proportional to the target pCO 2 conditions. Filtered seawater was introduced to each aquarium at a flow rate of 150 mL min −1 . The temperature of all experimental tanks was maintained at 17 • C and was slowly increased to 19.5 • C ( Supplementary Figure 1.2) on days 39-51 of the experiment in an effort to stimulate gonad development for a companion experiment. It should be noted that this is a small temperature shift relative to what these oysters experience seasonally in their native habitats ( Supplementary Figure 1.1). Oysters were fed 1% Shellfish Diet 1800 R twice daily following best practices outlined in Helm and Bourne (2004).
Temperature, pH, and salinity of all tanks were measured three times per week (M,W, and F) for the duration of the experiment. Seawater pH was measured with an Accumet solid state pH electrode (precision = 1 mV) calibrated with pH 7.01 and pH 10.01 NBS buffers (for calibration slope) and Dickson seawater Certified Reference Material (for calibration intercept). Complete carbonate chemistry was determined for each tank every 2 weeks. In brief, seawater samples were collected every 2 weeks in 250 ml borosilicate ground-glass stoppered bottles sealed with vacuum grease from each tank and immediately poisoned with 100 µl saturated HgCl 2 solution, then refrigerated until analysis of dissolved inorganic carbon (DIC) and total alkalinity (TA) was performed. DIC, TA, salinity, and temperature were used to calculate calcite saturation state, pH, CO 2− 3 , HCO − 3 , aqueous CO 2 , and pCO 2 of each sample using CO 2 SYS version 2.1 (Pierrot et al., 2011). Additional details about the experimental design and calculation of carbonate system parameters is provided in Supplementary Method 1.
Oyster EPF and tissue were sampled at six discrete time points (n total = 108, n = 6 per treatment per time point, see  : 24, 48 h, 9, 22, 50, and 80 days). This sampling schedule was designed to capture rapid physiological changes that occurred during the onset of exposure and to examine long-term stability of physiological and molecular responses.

Extrapallial Fluid Chemistry
Oyster pH EPF was measured by removing each oyster from its tank, inserting a 5 mL syringe with a flexible 18-gauge polypropylene tip through the luer-lock port into the oyster's extrapallial cavity, and extracting ∼0.5-2 mL of extrapallial fluid. Care was taken to avoid puncturing the mantle tissue and inadvertently sampling either the hemolymph or stomach fluid. The pH EPF was measured immediately after extraction with an Orion 91'10DJWP Double Junction micro-pH probe calibrated with pH 7.01 and 10.01 NBS buffers (for slope) and Dickson seawater Certified Reference Material (for intercept).

Calcification Rate
Net calcification rate was calculated for oysters surviving to either 50 or 80 days (n calcification = 11, 12, and 12 for control, 1,000 ppm, and 2,800 ppm treatments, respectively) from dry weights obtained at the start of the exposure and on day 33 or 34. Dry weights were estimated from buoyant weight measurements using the linear relationship between oyster dry weight and oyster buoyant weight derived empirically for oysters investigated in the present study (see Ries et al., 2009, for details). The average daily change in dry weight was divided by initial dry weight to standardize calcification rate for allometric effects, and multiplied by 100 to convert that fraction into a percent. See Supplementary Method 2 for a complete description of the buoyant weight measurements and the estimation of the dry weights.

Tissue Collection
Patterns observed in pH EPF over time informed the decision to focus molecular analyses on mantle edge tissue extracted on days 9 and 80 of the exposure in the control and high OA treatments (n tissue = 6 per treatment per time point). For consistency and to minimize potential effects on mantle tissue arising from drilling of the oysters' right valves when installing the EPF extraction ports, mantle tissue was only collected from the left valve of the oysters (i.e., the side opposite of the EPF extraction port). Immediately following extraction of EPF at each time point, the oysters were shucked and the edge of the mantle tissue was sampled. Tissue was immediately flash frozen in liquid nitrogen before being transferred to a −80 • C freezer for storage prior to DNA/RNA extraction.

DNA Methylation Library Preparation and Quantification
DNA was isolated from the oyster mantle edge tissue samples using the E.Z.N.A. Mollusc Kit (Omega) pursuant to the manufacturer's instructions. Isolated DNA was quantified, sonicated, and fragmentation was verified using a 2200 TapeStation System (Agilent Technologies, Santa Clara, CA, USA). Samples were enriched for methylated DNA using MethylMiner kit (Invitrogen, CA, USA), then sent to Zymo Research (Zymo,CA,USA) for bisulfite-conversion, library preparation, and sequencing. Libraries were sequenced on an Illumina HiSeq1500 platform with paired-end 100 bp sequencing.
Raw sequences were trimmed using TrimGalore! (v0.6.4; Martin, 2011) and aligned using Bismark (v0.22.3;Krueger and Andrews, 2011). Bismark was also used to determine the methylation status, methylated or unmethylated, for all sequenced CpGs. Files were then processed by performing a median based normalization of methylation counts among samples, destranding CpGs loci, and filtering loci that did not have at least 5x coverage for each sample using the R package methylKit (v1.10.0; Akalin et al., 2012). Finally, CpG loci were annotated by feature (e.g., introns, exons, and intergenic regions) and genomic coordinates using feature tracks created based on methods outlined in Venkataraman et al. (2020) and gene annotation file available on NCBI (Accession: GCA_002022765.4). Hereafter, "gene body" is used to reference the transcriptional region of the gene, which includes exons, introns, and untranslated regions, but not promoters. See the Supplementary Method 3 for a detailed description of this pipeline.

RNA Library Preparation and Quantification
RNA was extracted from the same mantle tissues used for DNA methylation using TRI Reagent Protocol (Applied Biosystems) following manufacturer's instructions. A cDNA library was constructed for each sample using a TruSeq Stranded mRNA kit (Illumina, Cat# RS-122-2101) following manufacturer's instructions. Library fragment size was determined using a BioAnalyzer (RNA 6000 Nano Kit; Agilent, CA, USA) and concentration was quantified on a Qubit (Qubit High Sensitivity RNA Kit; Invitrogen,CA,USA). Samples were sent to Genewiz (NJ,USA) for library preparation and sequencing. Libraries were sequenced on two lanes (12 individuals/lane) of an Illumina HiSeq 3000 platform with paired-end 200 bp sequencing.
RNA reads were trimmed for quality using Trimmomatic (v0.32; Bolger et al., 2014) and mapped to the oyster genome (Accession: GCA_002022765.4) using STAR (v.2.7.0; Dobin et al., 2013). Expression levels for each gene were quantified using RSEM (v1.3.2; Li and Dewey, 2011). Single sample gene count estimates were combined into count matrices (number of samples × number of genes) using a custom shell script (see Data Accessibility for additional script information), then filtered to remove genes that did not contain at least 1 count per million transcripts in at least 80% of individuals within at least one treatment-time experimental level. Expression data was then normalized using the "TMMwsp" method, transformed into log2 counts per million (log2-cpm), and sampled weighted in the R package "limma" (v3.40.6; Ritchie et al., 2015). See the Supplementary Method 4 for a more detailed description of this pipeline.

Extrapallial Fluid pH (pH EPF ) and Calcification Response
We tested the hypothesis that C. virginica pH EPF was associated with OA by using a linear mixed model. The full model included the explanatory variables of treatment, time, and their interaction as categorical fixed effects, and tank (nested within block) and oyster collection site as random effects. We conducted two analyses, one for each of the following response variables: (1) the measured value of pH EPF and (2) pH EPF relative to treatment pH seawater [ pH = pH EPF -pH seawater ; a metric used previously by Liu et al. (2020) as a better indicator of active pH EPF regulation]. Each model was performed in R using the lme4 package (v1.1-21; Bates et al., 2015). A stepdown strategy using likelihood ratio tests (based on the degrees of freedom estimated using the Satterthwaite method) was implemented in lmerTest, following the principle of marginality to select the most parsimonious model (v3.1-0; Kuznetsova et al., 2017). To evaluate significant differences in EPF chemistry (both pH EPF and pH) between the moderate and high OA treatments compared to the control treatment for each time point, post-hoc comparisons with a Tukey correction were performed using the multcomp package (v1.4-10; Hothorn et al., 2008). In addition, single sample t-tests were performed to test whether pH was significantly different from 0 (i.e., pH EPF was significantly different than pH seawater ) for oysters under each treatment at each time point. We corrected for testing multiple hypotheses using a Benjamini-Hochberg correction (Benjamini and Hochberg, 1995).
We examined the effect of pCO 2 on long-term EPF response (both measured pH EPF and pH) using linear mixed effect models. The full models included EPF response (pH EPF or pH) as the response variable, treatment and time and their interaction as categorical fixed effects, and tank (nested within block) and oyster collection site as random effects. Only samples from day-50 and day 80 time points were included in the model in order to investigate long-term trends in EPF response. Model selection and post-hoc comparisons were performed using the approach described in the previous paragraph.
We tested the hypothesis that calcification rate was associated with OA or pH EPF using linear-mixed effect models. The full models included calcification rate as the response variable, either pH EPF or treatment as a fixed effect, and tank (nested in block) and collection site as random effects. Time was not included as a fixed effect because calcification rate of oysters was measured across a single time interval (day 0 to day 33 or day 34). Treatment was handled as a continuous variable and was calculated as the average pCO 2 for each tank between days 0 and 33 or 34. Model selection was performed using the approach described in the first paragraph of this section. Regression analysis was used to examine the effect of either treatment or pH EPF on calcification rate.

Genome-Wide Molecular Response
To test the hypothesis that OA-induced changes in the transcriptome and methylome, we examined genome-wide gene expression and gene body DNA methylation using a principal components analysis (PCA) and PERMANOVA approach to visualize and test for differences amongst treatments and time points. The PERMANOVA was based on Manhattan distances and tested the null hypothesis of no effect of treatment, time, or their interaction on global gene expression and DNA methylation patterns using the R package vegan (v2.5-5; Dixon, 2003). To further investigate the stability of changes in genome-wide variation (gene expression or DNA methylation) induced by OA over time, we used a discriminant analysis of principal components (DAPC) performed with the R package adegenet (v2.1.2; Jombart et al., 2010). For our analysis DAPC was used to generate a discriminant function that maximally differentiated between treatment using samples from day 9, then we predicted where samples from day 80 would fall along this discriminant function to evaluate whether samples exhibited similar patterns of molecular differentiation across time. This allowed us to evaluate if changes that differentiated molecular variation among treatments at day 9 were similarly differentiated at day 80 (i.e., stable through time). See the Supplementary Method 5 for additional information on specific parameters used for the PERMANOVA and DAPC analyses.
Next, we further examined the hypothesis that OA induced changes in global DNA methylation by using a linear model to look at the effect of treatment, time, and feature (i.e., exon, intron, and intergenic region) on global methylation. Global methylation in this case was summarized as the median methylation across all CpGs within a feature for each sample.

Differential Molecular Response
The hypothesis that individual genes within the transcriptome were sensitive to OA was tested via differential gene expression analysis performed using a generalized linear model approach implemented in the R package limma (v3.40.6; Ritchie et al., 2015), using treatment, time, and their interaction as fixed effects. Site was not considered in this analysis given that it did not have a significant effect on either the phenotypic or genome-wide responses. Genes with FDR ≤ 0.05 and absolute value of log 2 fold ≥ 2 were considered differentially expressed.
We tested the hypothesis that OA induced changes in DNA methylation at individual CpGs (differentially methylated loci or DML) by performing a logistic regression using treatment as a fixed effect and time as a covariate implemented in the R package methylKit (v1.10.0; Akalin et al., 2012) using the default "slim" method to correct P-values. We also performed additional logistic regressions separately for each time point including treatment as a fixed effect. Only CpGs with coverage ≥5 for all samples were considered. Loci that had differential methylation of at least 50% between groups being compared and a q-value of <0.01 were considered significant. See the Supplementary Method 6 for additional information on the differential methylation and gene expression analyses.

Gene Co-expression Network Analysis
Clustering genes with similar levels of expression among individuals can help identify groups of genes that may be operating together, such as genes located within the same pathway. The response of these clusters to stress or stimulus can then be evaluated, which can potentially lead to the detection of transcriptomic responses that might otherwise be too subtle to be detected by standard gene-level differential expression analysis. We used a weight gene co-expression network analysis implemented in the R package WGCNA (v1.68; Langfelder and Horvath, 2008) to identify and cluster genes that co-express among individual oysters in order to reexamine the hypothesis that OA induced changes in gene expression and specifically evaluate whether particular clusters of genes respond to OA. We also used these gene clusters to evaluate the hypotheses that OA-induced changes in gene expression are (i) associated with shifts in DNA methylation and (ii) associated with regulation of the pH EPF . In brief, we followed a standard WGCNA pipeline which identified clusters of genes, termed modules, based on their level of dissimilarity among all individuals (regardless of treatment or time point) using the "Ward.D2" method for clustering (Murtagh and Legendre, 2014). Expression of each module for each individual was estimated as the eigenvalue (i.e., the first principal component value for each individual) from a PCA of an expression matrix generated from all genes included in the module for each individual. Similarly, DNA methylation was averaged for all CpGs located within genes included in a module to generate a single methylation value per individual. Linear regressions were used to identify gene expression modules that were significantly associated with EPF response (i.e., pH) and to test whether the expression of these modules was also associated with gene body methylation. See Supplementary Method 7 for additional details on the steps used in the WGCNA analysis and follow-up statistical analyses.

Functional Enrichment Analysis
As part of the hypothesis that OA induces changes in both gene expression and DNA methylation, we performed a functional enrichment test to investigate whether particular gene ontology (GO) categories were enriched in OA-induced molecular responses (either in gene expression or methylation) using GO-MWU, a rank-based gene enrichment method developed by Wright et al. (2015). We performed this analysis separately for each time point using the log 2 -fold change in gene expression and the difference in mean methylation among treatments. Mean methylation was calculated as the mean among all CpG loci within a gene across all individuals within a particular treatment and time point. Only genes with at least 5 CpG loci were considered for the analysis to ensure mean methylation estimates were based on genes with at least moderate CpGs coverage. We also used this analysis to determine if gene modules that were significantly associated with EPF (as determined in the previous section) were enriched in particular GO categories based on the module membership values for each gene within the module. These scores represent a measure of strength for each gene's membership to that module and are used here to determine if particular GO categories are enriched in each module. This analysis was run separately for GO categories associated with molecular function, biological process, and cellular components. A 10% FDR correction (GO-MWU default) was used to adjust for multiple comparisons. See Supplementary Method 8 for additional details on the GO-MWU analysis.

Gene-Level Characterization of Gene Expression and DNA Methylation
The hypothesis that DNA methylation correlates gene expression was tested using a PCA-based approach (described in Gavery and Roberts, 2013) to examine the relationship between gene expression, DNA methylation, and gene attributes (e.g., gene length and number of exons). The following variables were used in the input matrix for the PCA: mean gene expression over all treatments (mean log2-cpm values); the coefficient of variation (CV) in gene expression among treatment means; mean DNA methylation level over all treatments; the DNA methylation CV among treatment means, gene length, exon number per gene; and the number of CpGs per gene (scripts for generating each variable within the matrix are available on the github repository). Gene attributes were normalized by log transformation.

Extrapallial Fluid Chemistry and Calcification Response
We tested the hypothesis that C. virginica pH EPF was associated with OA using a linear model that included both fixed effects (treatment and time), their interaction, and the random effect of tank nested within block. There was a significant effect of treatment (F 2,14.783=5.8206 ; P = 0.0137) and the interaction between time and treatment (F 10,74.002=2.412 ; P = 0.0153) on pH EPF evaluated across all 6 time points. Post-hoc tests revealed a significant decrease in pH EPF within the high OA treatment compared to the control condition at the onset of the exposure (48 h), no difference on days 9-22, and a significant decrease on the final two time points (days 50 and 80; Figure 1A).
A slightly different result was found when we examined pH (pH EPF -pH seawater ) using a model that included both fixed effects (treatment and time), their interaction, and the random effect of tank nested within block. Both treatment (F 2,14.818=18.9487 ; P = 0.0001) and the interaction between treatment and time (F 10,74.036=2.0117 ; P = 0.0439) were significant, while the independent effect of time (F 5,74.050=0.935 ; P = 0.4637) was not significant (Figure 1B). Post-hoc comparisons between the control treatment and each OA treatment revealed a significant increase in pH under the high OA treatment relative to the control treatment at the onset of the exposure (9-22 days, Figure 1B, "Trt"). By days 50 and 80, pH was no longer statistically different amongst treatments, although mean pH in both OA treatments were higher than pH in the control. One-tailed t-tests showed that pH of oysters in all treatments was significantly lower than 0 (i.e., the pH EPF was more acidic than the pH seawater ) at almost all time points of the experiment (Figure 1B, "Env"), indicating a strong tendency for oysters to maintain a more acidic EPF fluid relative to their environment. The exception to the trend was the high OA treatment at day 9, where pH EPF was not significantly lower than pH seawater (Figure 1B).
We examined the long-term pH EPF and pH responses (i.e., conditions on days 50 and 80) using a model that included treatment as a fixed effect. Treatment was a significant predictor of both pH EPF (F 2,29 = 22.438, P < 0.0001) and pH (F 2,29 = 7.982, P = 0.0018; Figures 2A,B). The post-hoc tests revealed a significant decrease in pH EPF in the high OA treatment relative to both the control and moderate OA treatments, and showed an increase in pH in the two OA treatments relative to pH in the control.
Lastly, we found that calcification rate decreased in response to OA (R 2 = 0.409, P = <0.0001), which was driven largely by the substantial decline in calcification rate between the moderate and high OA treatments ( Figure 2C). The decline in calcification rate was also associated with a significant decrease in pH EPF (R 2 = 0.112, P = 0.028, Figure 2D).

DNA Methylation Responses to OA
Approximately 1.4 billion paired-end 80 bp reads of MBDenriched, bisulfite treated DNA were obtained across the all sequenced samples (NCBI BioProject ID: PRJNA594029, see Supplementary Table 3.1 for sequence coverage details). A total of 622.4 million quality filtered reads mapped to the C. virginica genome (88.3% of all CpGs in the genome, Supplementary Table 3.2). A single individual within the control treatment on day 9 was removed from downstream analysis due to poor sequencing depth and a high rate of gene duplication error after mapping (n tissue = 23). After filtering, we retained 403,976 CpGs with at least 5x coverage for each of the remaining 23 samples (Supplementary Table 3.2). Of the CpGs with at least 5x coverage, 94.1% were located within gene bodies, with 59.9% CpGs located in exons among 27,932 genes ( Figure 3A; "CpGs 5 "). CpGs with a minimum of 5x coverage in each sample were used for the downstream analyses.
We used a PERMANOVA to test the hypothesis that OA induced changes in genome-wide gene body methylation. We found significant differences in global methylation due to treatment (Adonis P Trt = 0.027), but not due to time (Adonis P Time = 0.235) or the interaction between treatment and time (Adonis P TimexTrt = 0.364; Figure 3B). This was supported by the DAPC, which found gene body methylation that was differentiated at day 9 was not maintained at day 80 ( Figure 3C). Moreover, global median methylation was significantly lower in the high OA treatment compared with the control treatment (hypomethylation, Figure 3D, P Trt < 0.0001), but significantly increased through time (hypermethylation, P Time = 0.001). However, this increase through time only led to a 1-5% shift in methylation when averaged over all loci in the genome. Feature  . OA", 1,000 µatm; "High OA", 2,800 µatm) was significantly different from the pH EPF in the control treatment ("Control", 580 µatm) for each time point (e.g., H 0 : Control pH EPF = High OA pH EPF ). (B) Symbols at the top of the graph indicate significant post hoc tests. The "vs. Ctrl" comparisons indicate time points where pH in one or both of the OA treatments was significantly different from the pH in the control treatment (e.g., H 0 : Control pH = High OA pH), while " pH = 0" comparisons indicate time points where pH in one of the treatments was significantly different from 0 (e.g., H 0 : Control pH = 0). The latter comparison is equivalent to evaluating whether pH EPF of a particular treatment is significantly different from its respective pH seawater . Statistical significance is denoted by asterisks ( *** P < 0.001, ** P < 0.01, * P < 0.05). Time points are presented on a non-linear scale and treatment points within time points were staggered along the x-axis to improve visualization.
Frontiers in Marine Science | www.frontiersin.org FIGURE 2 | Long term trends in pH EPF , pH, and calcification. (A) Long-term average pH EPF (total scale) by treatment; horizontal lines indicate mean seawater pH for the three treatments (dashed: control, dot-dashed: moderate OA, dotted: high OA). (B) Long-term average pH (pH EPF -pH seawater ; total scale) by treatment. (C) Percent daily calcification rate (%-change in shell mass/day) by treatment (slope = −1.751e-05, P < 0.0001, R 2 = 0.409), dotted line indicates calcification rate of zero percent and solid black line represents fitted regression. (D) Calcification rate vs. pH EPF (slope = 0.032, P = 0.028, R 2 = 0.112); solid black line represents fitted regression. pH EPF measured on day 50 or 80 (n = 35). Letters on first two panels (A,B) represent levels of significance based on post hoc testing, while bars on (C) represent standard errors.
had the greatest effect on global median methylation (P Feature < 0.0001), with exons exhibiting the highest levels of methylation.
The hypothesis that OA induced changes in the DNA methylation of specific CpGs was evaluated using a logistic regression approach implemented in methylKit. The initial model contained all samples and included treatment as the main effect and time as a covariate ( Figure 3A; ''DML Trt "). Only 2 CpGs were found to be differentially methylated by treatment, both located in gene bodies. Next, we looked at each time point separately including only treatment as the main effect ( Figure 3A; day 9-''DML Trt_D09 ;" day 80-''DML Trt_D80 "). This led to the additional discovery of 83 DMLs by treatment (25 DML on day 9 and 58 DML on day 80; Figure 3A). An additional 50 DML were significant through time in the control treatment, although there was little overlap between these DML and those responding to OA ( Figure 3E). Additionally, we saw no overlap of specific DML by treatment among days (Figure 3E), further illustrating that the loci induced by OA shifted through time as indicated by the DAPC (Figure 3C).
Of the 85 DML there was little functional overlap among DML at each time point, with the exception of several genes associated with the THO complex, a nuclear structure composed of multiple proteins involved in transcription elongation and mRNA maturation (Supplementary Table 3.3). On day 80, we found DML in several biomineralization and/or stress response related genes, including cadherin, protein ubiquitination, and death effector domain-containing genes. Two GO categories were enriched in hypomethylated CpGs, the cellular component "amidotransferase complex" on day 9 and the biological process "biosynthetic process" on day 80 (Table 1). No categories were shared among time points.

Gene Expression Responses to OA
RNA sequencing yielded a total of 955 million paired-end reads (NCBI BioProject ID: PRJNA594029) from all sequenced samples (see Supplementary Table 4.1 for read coverage and mapping details). Reads mapped to 37,098 genes, of which  Table 4.2). An exploratory PCA highlighted a single individual in the day 80 control treatment as an outlier compared to all other samples, which was subsequently removed from downstream analysis (n tissue = 23). The hypothesis that OA induced changes in genome-wide gene expression was evaluated with PERMANOVA. We found that treatment (R 2 = 0.051, Adonis P = 0.037) and time (R 2 = 0.067, Adonis P = < 0.001), but not their interaction (Adonis P = 0.214), had a subtle but statistically significant impact on genome-wide gene expression ( Figure 4A). The DAPC showed that the global patterns of gene expression that differentiated the two treatments at day 9 were not maintained at day 80, with the strongest directional change along the discriminant function in the day 80 samples occurring in the high OA treatment (Figure 4B).

20,387 remained after filtering (Supplementary
We also tested the prediction that the transcriptome in oysters is sensitive to OA by looking at the response of individual genes using differential expression analysis. Interestingly, no individual genes were significantly differentially expressed between treatments, including a number of genes known from the literature to be associated with biomineralization (Supplementary Table 4.3; Supplementary Figure 4.1). However, biomineralization-associated genes tended to be overexpressed relative to the mean expression level for all genes in both control and OA treatments ( Supplementary  Figure 4.2), which suggests that their activity is elevated within the mantle tissue. values were determined by using the DAPC package to estimate a function that maximally discriminates between treatment using samples from day 9, then the discriminant value for each sample from day 80 was predicted using the same discriminant function. Arrows indicate direction and degree of movement of expression patterns for control (blue) or OA (red) treatments along the discriminant function from day 9 to day 80. Both plots are based on log2-cpm gene expression.
We also tested whether OA induced changes in particular GO categories using a GO enrichment analysis implemented in GO-MWU. Several GO categories were found to be enriched in genes both up-and down-regulated within the high OA treatment relative to the control treatment ( Table 1), but few GO categories appeared to be shared among time points. The exceptions were two categories enriched in genes up-regulated ("DNA metabolic process" and "cellular nitrogen metabolic processes") and one enriched in genes down-regulated ("transcription factor complex") in the high OA treatment. Overall, more GO categories were enriched on day 80, including several metabolism and biosynthesis related biological processes. Several GO categories associated with oxidoreductase activity were also up-regulated on day 80.

Investigating Associations Between DNA Methylation, Gene Expression, and pH EPF
The following comparative analyses (next two sections) were performed using samples where gene expression and DNA methylation data overlapped (n tissue = 22). Only genes with at least 20% coverage of all CpGs were included in the analyses, which represents a reasonable trade-off between maximizing the number of retained genes in the analyses while removing genes with limited or no DNA methylation data (Supplementary Figure 3.1).

Associations Between DNA Methylation and Gene Expression
To test the hypothesis that DNA methylation is correlated with gene expression, we used a PCA that included measures of DNA methylation, gene expression, and other gene attributes. We found the first three PC axes explained >85% of the variance in the data (Figures 5A,B). PC1 was dominated by gene features, including a negative association with the length of the gene (PC axis loading, percent contribution to PC axis; gene_length: −0.504, 27%), the number of exons (exon: −0.454, 22%), the number of CpGs (totalCpG: −0.484, 25%). PC2 was loaded positively by mean methylation (−0.636, 40%), and negatively by methylation variation among treatments (Methylation_CV: 0.588, 33%), suggesting a negative relationship between percent methylation and the amount of variation in methylation among treatments. Finally, PC3 was dominated by gene expression (0.482, 23%), variation in expression among individuals (Gene_Expression_CV: −0.672, 46%), and methylation variation among treatments (Methylation_CV: 0.376, 15%), indicating a positive relationship between variation in methylation and gene expression and a negative relationship between gene expression and gene expression variation. Importantly, these patterns were robust even under increasingly stringent DNA methylation coverage scenarios (Supplementary Figures 3.2, 3.3). Finally, simple linear regression found that gene expression exhibited a significant positive relationship with gene body methylation (R 2 = 0.0624, P < 0.0001; Figure 5C), while the variation in expression had a significant negative relationship with gene body methylation (R 2 = 0.0658, P < 0.0001; Figure 5D).

OA-Induced Shifts in DNA Methylation, Gene Expression, and pH EPF
We tested the hypothesis that shifts in DNA methylation were associated with OA-induced changes in gene expression using regression analysis, with log2-fold change in gene expression as the response variable and the percent change in mean gene methylation (i.e., average DNA methylation for all CpGs with coverage in a gene) as a fixed effect. We found there was a significant positive effect of gene methylation changes, but the slope was very small and <2% of the variance in the log2-fold change in gene expression was explained by gene methylation changes (day 9-slope = 0.0121, P < 0.0001, R 2 < 0.012; day 80-slope = 0.0137, P < 0.0001, R 2 = 0.014; Figure 6A). This indicates that there was a weak but statistically significant relationship between DNA methylation and gene expression.
This weak relationship may have resulted from averaging methylation within each gene body, so regression analysis was used to examine the relationship between the change in methylation for each DML and the log2-fold change in expression of each corresponding gene. While we found a positive association between DML and corresponding log2-fold change in gene expression that explained a higher percent of the variation, the slope was substantially lower (slope Day9 = 0.0005, slope Day80 = 0.00054, P Day9 = 0.0005, P Day80 = 0.0004, R 2 Day9 = 0.406, R 2 Day80 = 0.276; Figure 6B) than in the previous gene-level methylation analysis. Thus, both analyses indicate that differential methylation had a substantially small effect on differential gene expression.
Gene co-expression network analysis identified 52 modules of co-expressed genes (i.e., genes that shared similar expression patterns among all 22 individuals) ranging from 36 to 5,059 genes in size (Supplementary Table 5.1). We used these modules to test the hypotheses that (i) pH EPF is associated with the expression of specific gene modules, (ii) the expression of these pH EPFassociated modules is also associated with DNA methylation, and (iii) that DNA methylation within these modules is sensitive to OA. Of the 52 modules, the eigengene expression of four modules was significantly associated with pH (excluding one module that was driven by a single individual outlier; Figure 7A). We did not see a significant association between the eigengene expression of any of these four modules and mean gene body DNA methylation (Figure 7B). Although, mean gene body DNA methylation within these modules was significantly associated with treatment in three of the four modules, but not time or the interaction of treatment and time ( Figure 7C). These modules were composed of a diverse collection of genes that were not significantly enriched in specific GO categories. There was, however, evidence of genes involved in pH regulation, ion transport, and other biomineralization-related processes included within the modules (Supplementary Table 5.1).

DISCUSSION
We conducted an integrative analysis of physiological, transcriptomic, and methylomic data in the mantle-edge tissue of eastern oysters exposed to a range of pCO 2 conditions to elucidate the oyster's capacity to regulate pH EPF in response to OA and examine the molecular underpinnings of this response. At the physiological level, oysters maintained a pH EPF significantly lower than pH seawater in our control treatments, which is consistent with what has been shown previously in oysters and other species of bivalves (e.g., Crenshaw, 1972;Cameron et al., 2019). This trend was also seen in the OA treatments at most time points, the exception being exposure day 9 in the high OA treatment where oysters maintained pH EPF near pH seawater . On day 9 and 22 of the exposure, oysters in the high OA treatment also maintained a pH significantly higher than pH in the control, which resulted in a non-significant difference in the pH EPF between the high OA and control treatment on those days.
At a molecular level, C. virginica exhibited a subtle genomewide response in both the methylome and transcriptome to OA, but surprisingly little response in individual gene expression or methylation. When genes were clustered, we found evidence of a significant association between gene clusters and pH, suggesting subtle shifts in these genes could be responsible for regulating pH EPF in response to OA. However, gene body methylation within these gene clusters was not associated with gene expression, indicating that DNA methylation did not have a significant role in gene regulation. Below, we evaluate the short and long-term effects of OA on both phenotypic and molecular response, examine how these results link together and relate to the current model of EPF fluid regulation, and discuss possible explanations for the surprisingly subtle shifts in both the transcriptome and methylome.

pH EPF Response Through Time
Crassostrea virginica, like other coastal calcifiers, are experiencing OA due to a combination of local (Wallace et al., 2014) and global change-associated processes (e.g., Orr et al., 2005). In particular, local processes, such as excessive nutrient loading (Wallace et al., 2014), input from acidic river water (Salisbury et al., 2008), and upwelling of CO 2 rich groundwater (Feely et al., 2008), are exposing oysters to more intense short-term acidification events well before global, anthropogenic CO 2 -induced acidification is predicted to impact these systems. Consistent with these findings, water chemistry data from our collection sites found substantial variation in pCO 2 , including intervals of acidification FIGURE 6 | Gene expression change in response to OA as a function of DNA methylation change. Log2-fold change in gene expression compared to % difference in (A) mean gene methylation among treatments (Day 9 -slope = 0.0121, P < 0.0001, R 2 = 0.012; Day 80 -slope = 0.0137, P < 0.0001, R 2 = 0.014) and (B) significantly differentially methylated loci (DML) among treatments (Day 9 -slope = 0.0005, P = 0.0005, R 2 = 0.406; Day 80 -slope = 0.00054, P = 0.0004, R 2 = 0.276). Orange line represents the fitted linear model. Colors in (B) correspond to DMLs that were significantly hyper-(green) and hypomethylated (purple) in the OA treatment.
near the high OA scenario employed in the present study. Despite the apparent frequency of these short-term acidification events, few studies have examined the response of adult oysters to such high levels of OA (i.e. the high OA treatment in present study), relative to the lower levels of OA-associated with increasing anthropogenic pCO 2 (i.e. the moderate OA treatment in present study).
Bivalve response to OA has received considerable attention over the past two decades (see review Gazeau et al., 2013), with a substantial focus on measuring calcification (Miller et al., 2009;Talmage and Gobler, 2009;Waldbusser et al., 2011a) and physiological responses (Beniash et al., 2010;Ivanina et al., 2013;Matoo et al., 2013;Gobler and Talmage, 2014) in C. virginica. Recent research has shown that bivalves can isolate and moderate the fluid at the site of calcification, and elevate the pH of this fluid relative to seawater under OA conditions (Cameron et al., 2019;Liu et al., 2020). It has been hypothesized that this may be a strategy for maintaining calcification even under increasingly acidic conditions. However, the degree to which a species can moderate its calcifying fluid pH under OA scenarios varies across taxa (e.g., McCulloch et al., 2012;Holcomb et al., 2014;Liu et al., 2020), life stage (e.g., Ramesh et al., 2017), and temperature (e.g., Cameron et al., 2019). Previous work with C. virginica found that they exhibited a moderate ability to regulate their pH EPF after an extended 60-day experimental exposure, elevating their pH EPF above pH seawater under most OA scenarios (Liu et al., 2020).
We predicted that similar to this previous study, C. virginica would increase its pH EPF relative to pH seawater in response to OA. Consistent with our prediction and with previous studies on C. virginica (Liu et al., 2020) and adult King scallops Pectens FIGURE 7 | Co-expression and DNA methylation responses to OA. Four of the top five modules associated with oyster pH (colors are randomly assigned by the WGCNA package and have no meaning). Scatterplots of (A) pH by eigengene expression (P cyan = 0.021, P lavenderblush3 = 0.044, P darkred = 0.037, P turquoise = 0.042), (B) eigengene expression by mean module methylation (P cyan = 0.087, P lavenderblush3 = 0.052, P darkred = 0.264, P turquoise = 0.833), and boxplot of (C) mean module methylation by treatment and time point. For (C) there was a significant effect of treatment in three of the four modules (P cyan = 0.035, P lavenderblush3 = 0.038, P darkred = 0.123, P turquoise = 0.013), but neither time point or their interaction was significant for any of the modules. Mean module methylation was calculated as the mean methylation for all CpGs for each gene within a module. Solid lines indicate a significant relationship between explanatory (x-axis) and response variable (y-axis) using a linear model, while dotted lines indicate non-significant trends among variables.
maximus (Cameron et al., 2019), we observed an increase in pH under OA conditions. However, pH EPF remained below pH seawater at most time points under the different treatments, which is inconsistent with the prior work on C. virginica (Liu et al., 2020) and P. maximus (Cameron et al., 2019). One possible explanation for the lower pH EPF and pH in the present study compared with Liu et al. (2020) was the choice of experimental temperatures (∼18-19 • C in the present study vs. 25 • C in the other study), which has been shown in P. maximus to impact the relationship between pCO 2 and pH EPF (Cameron et al., 2019). The lower pH EPF in the present study may also reflect increased respiration of CO 2 into the EPF, potentially driven by increased metabolic activity and/or rates of feeding (Crenshaw, 1972;Sutton et al., 2018;Cameron et al., 2019). Comparing among treatments, we found that pH was significantly higher in our high OA treatments relative to the control at days 9 and 22, indicating that OA prompted a response in pH EPF . This shortterm response could be the result of stress-induced metabolic depression or increased removal of hydrogen ions from the EPF via proton transporters or other ion regulators. Moreover, the short-term elevation of pH in OA treatments was not sufficient to prevent net shell loss (i.e., net dissolution) after 35 days of exposure. However, a prior study on C. virginica suggests that gross calcification does occur under comparably high OA conditions, despite net shell loss . Future studies that specifically consider the deposition rate of calcium carbonate and the impact of the interaction between temperature and OA on the regulation of EPF chemistry in C. virginica and other bivalves will offer substantial insight into the importance of EPF regulation on this species' calcification response to OA.

Biomineralization Response
Recent studies have begun to clarify the role of particular genes, pathways, and associated molecular mechanisms in the biomineralization process, and their responsiveness to OA across a diverse group of marine calcifiers (see review by Clark, 2020). In bivalves, a principal component of the biomineralization process is the transport of ions between the mantle epithelium and the EPF, in particular Ca 2+ and HCO − 3 , needed for the formation of calcium carbonate. The movement of Ca 2+ ions into the EPF occurs both passively via paracellular pathways and by active transport through transporters in the mantle epithelium membrane (primarily Ca 2+ -ATPase and Na + /Ca 2+ exchangers), while HCO − 3 is shuttled into the EPF via HCO − 3 transporters. Additionally, carbonic anhydrase is an important protein in the biomineralization process that catalyzes a reaction that generates HCO − 3 and H + ions from CO 2 and water. The H + ions produced in this reaction (among many others) serve a critical role in several of the transporter proteins both directly (e.g., Ca 2+ -ATPase) and indirectly (e.g., Na + /H + exchanger) involved in biomineralization. In addition, reducing the amount of H + ions in the EPF (measured as the pH of the fluid) causes HCO − 3 to dissociate to CO 2+ 3 , thereby increasing the calcium carbonate saturation state in support of calcification (Ries, 2011). A number of other marine calcifiers have also exhibited shifts in genes associated with these pathways of ion-transport in response to CO 2 -induced OA (Evans et al., 2013;Davies et al., 2016;Johnson and Hofmann, 2017;Griffiths et al., 2019).
We hypothesized that OA would elicit a response in genes involved in biomineralization in the mantle, particularly those involved in ion-transport. Based on the shift in pH from days 9 to 80, we predicted a similar shift in biomineralizationrelated genes between these same time points. However, we found little evidence of differential gene expression or methylation corresponding to biomineralization or ion transport between these days or across treatments. Despite the lack of differentiation among treatments, biomineralization genes were well-represented in our data and overexpressed relative to all other genes, consistent with current biomineralization models in bivalves. Although carbonic anhydrase and several ion transporter genes were not significantly differentially expressed, they did display a subtle up-regulation in response to OA on day 9 relative to day 80 (Supplementary Figure 4.1), indicating that these subtle responses in biomineralization genes could have contributed to the significant global response to treatment that was observed. These results are surprising given the more pronounced gene expression responses seen in other marine calcifiers (Evans et al., 2013;Davies et al., 2016;Johnson and Hofmann, 2017;Griffiths et al., 2019). One possible explanation for the surprisingly weak response observed in the present study was the decision to examine a single tissue type, the mantle, rather than look at the transcriptomic response of the entire organism or multiple tissues (but see Hüning et al., 2013;Li et al., 2016). Moreover, these subtle responses appear to be fairly consistent with previous targeted gene expression studies examining C. virginica mantle tissue in response to OA exposure, which also did not find significant differential expression of biomineralization genes (Beniash et al., 2010;Richards, 2018). Thus, biomineralization-associated pathways in C. virginica mantle tissue may be relatively canalized compared to other marine calcifiers.

Stress and Metabolic Responses
A broad range of transcriptomic responses related to stress and metabolism have been noted in marine calcifiers in response to OA (see review by Strader et al., 2020). In particular, OA exposure has been shown to elicit an up-regulation of oxidoreductaserelated genes and oxygen-stress related enzymes in a broad range of organisms, including corals (Kaniewska et al., 2012), sea urchins (Todgham and Hofmann, 2009), and oysters (Tomanek et al., 2011;Wang et al., 2016). In C. gigas, OA exposure led to a delayed antioxidant and oxidative stress response in gill tissue, suggesting that redox homeostasis was possible for short-term (up to 14 days) but not long-term (28 days) exposure to OA (Wang et al., 2016). Consistent with this response in C. gigas and paralleling the shift we observed in the partial maintenance of the pH EPF , we found several GO categories associated with metabolic response, in particular oxidative stress. In general, we observed a larger stress response over time, suggesting that OA induced a mild stress response in the mantle tissue, which was exacerbated over time. We did not find evidence of overlap between these GO categories in gene expression and DNA methylation, but we did observe DML within genes corresponding to a number of different stress responses, including hypomethylation of death-effector domain-containing genes and hypermethylation of a protein-ubiquitination related gene on day 80. Both types of genes are involved in cell apoptosis and possibly reflect a stress response to OA. Similar OA-induced responses in protein-ubiquitination genes were found in the gonadal tissue of C. virginica (Venkataraman et al., 2020), suggesting OA can induce parallel shifts in DNA methylation in the methylome of different tissues. However, general overlap of OA-induced DNA methylation was low among tissues (2 OA-induced DMLs among 85 identified in the mantle were also induced in the gonadal tissue; Supplementary Figure 3.4). Future studies focusing on OA-induced changes in gene expression and DNA methylation among tissues may help clarify if the subtle responses in C. virginica are the result of a relatively canalized response in the mantle tissue or represents a pattern unique to the species.
The Molecular Basis for pH EPF Response to OA Ries (2011) proposed a model to describe the regulation of the calcifying fluid whereby the active pumping of protons out of the calcifying fluid elevates the calcifying fluid pH and maintains a favorable calcium carbonate saturation state ( > 1) even under highly acidic conditions. Empirical support for this model within oysters has been provided from shell-based boron-isotope estimates of oyster EPF (e.g., Sutton et al., 2018;Liu et al., 2020), but little is known about the molecular basis of this response. To test the predictions proposed by this model, we hypothesized that regulation of oyster pH EPF would be driven by changes in gene expression, particularly genes associated with biomineralization, ion transport, and proton pumping. In addition, we also hypothesized that OA-induced shifts in the transcriptome of C. virginica would correspond to shifts in the methylome. This was based on growing evidence in the DNA methylation literature indicating that DNA methylation is linked with gene regulation (Zilberman et al., 2007;Zemach et al., 2010;Gavery and Roberts, 2013;Olson and Roberts, 2014;Dixon et al., 2018;Zhang et al., 2020), is induced by OA and other environmental stressors (Putnam et al., 2016;Dixon et al., 2018;Bogan et al., 2020), and corresponds with calcification response to OA in marine invertebrates (Liew et al., 2018).
By clustering genes based on their co-expression, we found that gene expression was associated with EPF response (i.e., pH) in four separate gene modules. This result, coupled with the lack of differential expression in individual genes, suggests that regulation of EPF occurs through subtle changes in a large number of genes. Interestingly, these correlated clusters were composed of a diverse collection of genes which were not significantly enriched in specific GO categories, although they did include genes involved in pH regulation, ion transport, and other biomineralization related processes.
Consistent with previous findings in marine invertebrates (Olson and Roberts, 2014;Liew et al., 2018), gene body DNA methylation in C. virginica mantle tissue exhibited a subtle but significant association with both gene expression (positive) and gene expression variation (negative). However, the small R 2values and slopes from the regressions indicate that changes in gene body DNA methylation had only a small effect on changes in gene regulation in C. virginica mantle tissue. Moreover, we found there was a weak positive association between OAinduced log2-fold changes in gene expression and shifts in gene body DNA methylation. Although, we found there was no overlap among enriched GO categories based on gene expression compared to DNA methylation data. We also found no evidence of an association between the eigen gene expression of the four pH-associated modules and gene body DNA methylation, which suggests that an alternative mechanism might be regulating the expression of genes within those modules. Follow-up investigations that use a more targeted, gene-specific approach would clarify whether these weak signals represent an important part of oyster EPF regulation in response to OA.

Subtle Shifts and Constraints
The subtle transcriptomic response of C. virginica found in this study contrasted the more pronounced transcriptomic responses exhibited by other marine calcifying taxa in response to OA (e.g., Li et al., 2016). However, other examples of subtle gene expression shifts in response to environmental change do exist (Malik et al., 2020), and these minor shifts can still have phenotypic consequences (Wu et al., 2010). In C. virginica, these subtle molecular responses may be the result of: (i) extended evolution under regular extreme OA conditions; (ii) constrained plasticity of the mantle (e.g., an evolved canalization of gene expression in response to OA); and/or (iii) a highly polygenic genetic response architecture.
Crassostrea virginica inhabit coastal estuaries, which can experience regular, unpredictable, and rapid swings in acidity leading to conditions near or above the high OA scenario (Waldbusser and Salisbury, 2014). Since the high OA treatment is not unusual from the oyster's perspective, the subtle molecular responses we observed may reflect that the prescribed high OA treatment in the present study was not different enough from what oysters typically experience on diurnal and/or seasonal cycles to elicit a strong molecular response. On the other hand, highly uncertain environments, in which the frequency of fluctuations is much higher than what the organism can predict or respond to, can lead to the evolution of a single expressed phenotype (in this case, gene expression in the mantle) across all environments that maximizes the geometric mean fitness (Seger and Brockmann, 1987;Starrfelt and Kokko, 2012;Botero et al., 2015). The variable environments inhabited by the oysters in the present study may consequently select for constrained (e.g., canalized) responses in the transcriptome and methylome, in which the expression profile results from the evolution of a nearly fixed bet-hedging strategy, rendering the genome capable of only subtle shifts. From an energetic cost perspective, one would expect selective pressure against increases in mRNA and protein expression due to the costs associated with increased expression (Wagner, 2005;Weiße et al., 2015). Interestingly, calcification genes were associated with higher expression than the genome-wide background in the mantle tissue, suggesting that their expression may be maximized, which could constrain their ability to be further up-regulated in response to OA. Lastly, while we cannot completely exclude the possibility that prior acidification events in the field acclimated the oysters to OA prior to the start of the experiment, thereby limiting differences among treatments. This seems unlikely to completely explain our lack of differences since we performed a one month acclimation at control conditions prior to the experimental exposure. Moreover, if it is true that field-collected oysters show subtle, but significant physiological responses, then the studies that focus on organisms reared in control conditions may be biased toward larger responses than actually exist in nature.
Subtle shifts may also be observed if the traits responding to OA are highly polygenic (Pritchard et al., 2010) and/or there is a redundant genetic architecture (Yeaman, 2015). Although this body of literature is focused primarily on shifts in allele frequencies in response to selection rather than gene expression, the same principles about polygenic traits apply. In the case of plastic polygenic traits, the magnitude of response of individual genes may be relatively small and, therefore, not detectable via differential expression analysis. In the case of redundant genetic architecture, theory predicts that high genetic redundancy in a population occurs when individuals can reach a phenotypic outcome by a wide range of genes and pathways (see topic reviewed in Láruson et al., 2020), which could result in very large residual error and limited statistical power to detect differences amongst treatments (i.e., lack of a clear significant response because every individual does something different to achieve the same phenotype). Regardless of whether these results are more influenced by subtle shifts or constraints on gene expression, they show that a marine calcifier's drastic phenotypic responses to OA (i.e., pH EPF response and reduced calcification rate under high OA) can be associated with only subtle molecular responses within that organism's calcifying tissue.

CONCLUSION
Crassostrea virginica exhibited an overall negative response to long-term OA (e.g., reductions in calcification and pH EPF accompanied by the subtle emergence of stress related transcriptomic responses). Although elevation of pH under the high OA treatment was enhanced at day 9, pH EPF was highly variable within all treatments during the first 22 days of the experiment-which indicates C. virginica may have a limited capacity to continue calcifying during ephemeral acidification events, such as those they already experience due to diurnal and seasonal fluctuations in seawater pH. Notably, we did not find a strong association between specific DML and differentially expressed genes or pH, suggesting that the regulation of pH EPF may occur by more subtle or diverse shifts in the methylome and transcriptome or via alternative mechanisms. Since the majority of OA transcriptomic studies have not specifically targeted calcifying tissue, our results highlight a need for more studies to focus specifically on calcifying tissues in order to assess how flexible the bivalve calcification process is at the molecular level, particularly in oysters. In contrast to other studies, we did not find strong evidence that DNA methylation acted as a mediator to gene expression or EPF responses to OA, suggesting that the importance of DNA methylation may vary considerably among tissues within the same species, among studies, and among taxa. If the subtle shifts we observed in this study are common in the calcifying tissue of other taxa, it may be difficult to gain an unbiased understanding of this phenomenon due to publication bias against weak or non-significant results. Thus, our study provides a roadmap for others seeking to understand how the calcifying fluid is regulated at the molecular level, even when patterns in the resulting molecular data are relatively weak. Subtle shifts at the molecular level and the capacity for marine calcifiers to maintain these molecular changes provide important context for understanding the capacity of different taxa to respond to long-term OA.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI BioProject ID -PRJNA594029. Scripts and non-sequence data can be found at the online github repository: https://github.com/ epigeneticstoocean/AE17_Cvirginica_MolecularResponse. The data can also be found on Dryad with link: https://doi.org/10. 5061/dryad.8cz8w9gnk.

AUTHOR CONTRIBUTIONS
Principle design of the experiment was done by AD-W, SR, JR, and KL. Field collections were done by AD-W, BF, and EM. The experiment was conducted by AD-W. During the experiment, extrapallial fluid was collected and analyzed by LC, while other phenotypic data and tissues were collected by AD-W. Water chemistry was analyzed by AD-W and EM. Tissues were prepared for RNAseq sequencing by BF and prepared for MBD-BSseq by YV and SR. All bioinformatic steps and statistical analyses were performed by AD-W. The manuscript was written by AD-W with assistance from KL, LC, JR, and YV. All authors contributed to editing and approving the final version of the manuscript.

FUNDING
This research was funded by the National Science Foundation (1635423). Resources purchased with funds from the NSF FMSL program (DBI 1722553, to Northeastern University) were used to generate data for this manuscript.

ACKNOWLEDGMENTS
We acknowledge Sara Schaal, Aki Laruson, Isaac Westfield, and other members of the Lotterhos and Ries labs for their helpful comments that improved the manuscript. This research was funded by the National Science Foundation (BIO-OCE 1635423). Resources purchased with funds from the NSF FMSL program (DBI 1722553, to Northeastern University) were used to generate data for this manuscript. We are also thankful for the help of many graduate students, undergraduates, technicians, and volunteers during the experiment, including Andrea Unzueta Martinez, Bodie Weedop, Isabel Gutowski, Jaxine Wolfe, Sandi Scripa, Kevin Wasko, Isaac Westfield, and Camila Cortina.