Changes in Genome-Wide Methylation and Gene Expression in Response to Future pCO2 Extremes in the Antarctic Pteropod Limacina helicina antarctica

Epigenetic processes such as variation in DNA methylation may promote phenotypic plasticity and the rapid acclimatization of species to environmental change. The extent to which an organism can mount an epigenetic response to current and future climate extremes may influence its capacity to acclimatize or adapt to global change on ecological rather than evolutionary time scales. The thecosome pteropod Limacina helicina antarctica is an abundant macrozooplankton endemic to the Southern Ocean and is considered a bellwether of ocean acidification as it is highly sensitive to variation in carbonate chemistry. In this study, we quantified variation in DNA methylation and gene expression over time across different ocean acidification regimes. We exposed L. helicina antarctica to pCO2 levels mimicking present-day norms in the coastal Southern Ocean of 255 μatm pCO2, present-day extremes of 530 μatm pCO2, and projected extremes of 918 μatm pCO2 for up to 7 days before measuring global DNA methylation and sequencing transcriptomes in animals from each treatment across time. L. helicina antarctica significantly reduced DNA methylation by 29–56% after 1 day of exposure to 918 μatm pCO2 before DNA methylation returned to control levels after 6 days. In addition, L. helicina antarctica exposed to 918 μatm pCO2 exhibited drastically more differential expression compared to cultures replicating present-day pCO2 extremes. Differentially expressed transcripts were predominantly downregulated. Furthermore, downregulated genes were enriched with signatures of gene body methylation. These findings support the potential role of DNA methylation in regulating transcriptomic responses by L. helicina antarctica to future ocean acidification and in situ variation in pCO2 experienced seasonally or during vertical migration. More broadly, L. helicina antarctica was capable of mounting a substantial epigenetic response to ocean acidification despite little evidence of metabolic compensation or recovery of the cellular stress response in this species at future pCO2 levels.


INTRODUCTION
Marine ecosystems are already reaching extremes of environmental change on par with projected global climate change over the next century (Hoegh-Guldberg and Bruno, 2010;Harris et al., 2013;Chan et al., 2017;Oliver et al., 2018). Species inhabiting regions such as polar oceans are currently faced with a necessity to acclimatize or adapt to stressors that will only intensify as climate change progresses. Driven by the rapid advancement of extreme physical conditions in ecosystems today and in modeled projections, global change biologists have begun to direct attention to the mechanisms and consequences of species' abilities to rapidly respond to anthropogenic stress via acquired, adaptive traits (Chown et al., 2007;Calosi et al., 2008;Chevin et al., 2010;Nicotra et al., 2010;Beldade et al., 2011;Donelson et al., 2018;Kelly, 2019) and how mechanisms of acclimatization vary within a system between current and projected environmental change (Hennige et al., 2010;Duarte et al., 2018). We examined such a rapidly-acting processnamely, changes in DNA methylation and gene expression in the Antarctic pteropod Limacina helicina antarctica during exposure to present and future ocean acidification (OA) levels for the Southern Ocean.
In recent studies of marine metazoans, epigenetic processes have been demonstrated to be associated with the phenotypic plasticity of fitness-related traits in species experiencing climate change-driven stressors (Zhang et al., 2013;Putnam et al., 2016;Clark et al., 2018;Strader et al., 2019;Wong et al., 2019), a relationship that can be explained in part by the link between epigenetic mechanisms, gene expression, translation, and the traits underpinned by these processes (True et al., 2004). For example, experimental manipulation of DNA methylation levels in Arabidopsis thaliana has been shown to substantially alter the phenotypic plasticity of key developmental and physiological traits in low-and high-nutrient environments (Bossdorf et al., 2010). In humans, quantitative trait loci at methylated bases are linked to differential gene expression underpinning multiple physiological traits across a diversity of tissue types (Taylor et al., 2019). Overall, understanding how epigenetic processes may promote adaptive responses to environmental stressors will require investigations in diverse, ecologically critical taxa and experimentation across gradients of stress eliciting adaptive and pathological responses.
Among a suite of epigenetic modifications, DNA methylation, the addition of a -CH 3 methyl group to either cytosine or adenine bases, has received considerable focus for its role in regulating gene expression, particularly in vertebrate lineages and plants (Schubeler, 2015;Zhang et al., 2018). Investigations into the role of DNA methylation in invertebrates have uncovered marked differences in methylation's influence on gene expression and performance between phylogenetic lineages and across different environmental contexts (Sarda et al., 2012;Dimond and Roberts, 2016;Hofmann, 2017;Eirin-Lopez and Putnam, 2019). However, studies of DNA methylation and its role in environmentorganism interactions in marine invertebrates remain sparse and confined to a small number of species for any given phyla.
Antarctic pteropods offer a valuable system for examining the role of epigenetics in shaping organismal responses to global change in the marine environment. Shelled pteropods or thecosomes have been proposed as a bellwether species for OA (Manno et al., 2017). The Antarctic thecosome L. helicina antarctica is a widely distributed and abundant pteropod endemic to the Southern Ocean (Boysen-Ennen et al., 1991;Steinberg et al., 2015;Thibodeau et al., 2019) where it can make up > 50% of total zooplankton measured as individuals per unit volume (Hunt et al., 2008). Limacina sp. is also a globally distributed genus spanning both polar and temperate environments (Bernard and Froneman, 2009;Hunt et al., 2010). L. helicina antarctica presently experiences a large degree of seasonal variability in seawater pCO 2 and pH, with pH varying by 0.6 units in regions like Prydz Bay (Gibson and Trull, 1999) and the Ross Sea (McNeil et al., 2010). This seasonal variability currently exceeds modeled predictions of an increase in annual mean pH of 0.4 units within the Southern Ocean by 2100 (IPCC, 2013). L. helicina antarctica also experiences present-day undersaturation in aragonite in some regions of the Southern Ocean during the austral winter (Bednaršek et al., 2012b). The frequency and duration of undersaturation events are expected to increase as OA advances in the Southern Ocean (McNeil and Matear, 2008;Hauri et al., 2015;Negrete-García et al., 2019). The onset of month long undersaturation events are predicted in the Southern Ocean within 14-29 years, along with annual mean aragonite undersaturation as soon as 10-20 years from today (IPCC, 2019).
In the face of future increases in pCO 2 , current evidence suggests that L. helicina antarctica is poorly poised to acclimatize to predicted OA over the next century (Bednaršek et al., 2012a(Bednaršek et al., , 2014bGardner et al., 2018). Juvenile L. helicina antarctica collected from the McMurdo Sound, Ross Sea, do not appear to alter their metabolic rate over the course of 2-week acclimations to a regional, present-day pCO 2 extreme of 427-513 µatm (Hoshijima et al., 2017). By contrast, acclimation to future pCO 2 extremes of 901-1000 µatm has been shown to significantly alter metabolic rate under near-ambient temperatures (Seibel et al., 2012;Hoshijima et al., 2017). Similarly, gene expression by L. helicina antarctica has shown little differential expression in response to intermediate pCO 2 levels of 432 µatm when compared to differential expression under 902 µatm at which pervasive downregulation is apparent among transcripts even after 21 days of acclimation (Johnson and Hofmann, 2017). The dynamic shifts in metabolism and gene expression that arise when L. helicina antarctica are exposed to present-day vs. future pCO 2 extremes may be driven or regulated by epigenetic mechanisms such as DNA methylation.
Determining (i) whether DNA methylation in L. helicina antarctica varies across exposures to current and future pCO 2 extremes and (ii) whether changes in DNA methylation are associated with differential gene expression may provide valuable insight into molecular mechanisms underpinning this species' performance under OA and its ability to rapidly acclimatize to changing ocean conditions. Findings regarding epigenetic responses by L. helicina antarctica to OA may provide translational information for similar studies in other members of Limacina sp. and thecosome pteropods, a clade of shelled pteropods considered to be a bellwether for the severity of global OA (Bednaršek et al., 2014a). Furthermore, epigenetic studies revealing the importance of DNA methylation for biological functions in molluscs (Gavery and Roberts, 2013;Riviere et al., 2013Riviere et al., , 2017Diaz-Freije et al., 2014;Lian et al., 2015;Garcia-Fernandez et al., 2017;Suarez-Ulloa et al., 2019) thus far have been heavily skewed toward bivalves and in particular, oysters, necessitating a broader assessment of diversity in environmental epigenetics across the speciose molluscan phylum.
To this end, we conditioned juvenile L. helicina antarctica collected from McMurdo Sound under present-day pCO 2 levels of 255 µatm, present-day pCO 2 extremes of 530 µatm, and future pCO 2 extremes of 918 µatm, resulting in aragnoite undersaturation, for up to 7 days (Figure 1). Following this experiment, we quantified the proportion of 5-methylcytosine in L. helicina limacina DNA following 1, 3, and 6 days of exposure to each pCO 2 level. Additionally, we sequenced and analyzed the transcriptomes of pteropods sampled from each pCO 2 level after 0.5 and 7 days of conditioning in order to (i) quantify differential expression in response to current and future pCO 2 extremes across time and (ii) test for correlations between differential expression and signatures of DNA methylation among transcribed genes.

Collection and Conditioning of L. helicina antarctica
Juvenile L. helicina antarctica were collected in October of 2014 from Cape Evans, McMurdo Sound, via plankton net as previously described (Hoshijima et al., 2017;Johnson and Hofmann, 2017) and transported to the Crary Laboratory, McMurdo Station, where they were held under near-ambient temperature (−0.59 • C) in filtered seawater for 24 h before conditioning at different pCO 2 levels. During conditioning, L. helicina antarctica were held in 1 L tanks with a flow rate of 2 L h −1 at a starting density of 200 animals L −1 . L. helicina antarctica were not fed at any point during the experiment as the wintersummer transition in the McMurdo Sound is characterized by low phytoplankton abundance (Foster et al., 1987).
Seawater pCO 2 was manipulated using three reservoir containers held at different pCO 2 levels, each connected to three replicate culturing tanks of the volume and flow rate described above. pCO 2 was controlled inside of reservoir tanks by injecting mixtures of filtered, CO 2 -scrubbed, dry air and pure CO 2 with seawater. Pure air and CO 2 were mixed using SmartTrak R 100 Series Mass Flow Controllers and Micro-Trak © 101 Series Mass Flow Controllers (Sierra Instruments, United States), respectively. Reservoirs were held in a 1240 L seawater table filled with near-ambient seawater. Culture tanks held a mean temperature of −0.59 • C ± 0.11 SD over the course of the 7-day experiment. Seawater temperature, pCO 2 , pH, and carbonate chemistry were measured as previously described (Johnson and Hofmann, 2017). Low, intermediate, and high pCO 2 levels averaged 255 µatm ± 3.52 SD, 530 µatm ± 12.03 SD, and 918 µatm ± 19.53 SD. A full report of temperature and seawater chemistry in each culture tank over time is available in Supplementary Figure S1.
Quantification and Analysis of Genomic 5-Methylcytosine treatment at 1, 3, and 6 days using a CTAB DNA extraction modified from Worden (2009). The step-by-step DNA extraction protocol is detailed under Supplementary Data S2. DNA purity, quality, and concentration were respectively assessed via Nanodrop, agarose gel, and Qubit Broad Range DNA quantification. DNA extractions were performed on five biological replicates per pCO 2 level x timepoint by sampling two pooled replicates from two of the triplicate cultures and one pooled replicate from the third.
The proportion of 5-methylcytosine (5-mC) present in L. helicina antarctica DNA was quantified colorimetrically using the ELISA-based MethylFlash DNA 5-mC Quantification Kit (Epigentek, Farmingdale, NY, United States) according to manufacturer guidelines. Contrary to manufacturer instructions, standards and samples were measured in triplicate rather than duplicate in order to improve accuracy. Additionally, the 45 samples were randomly split between two different plates such that 3/5 or 2/5 of replicates for each pCO 2 level x timepoint group would be present on a given plate. Within a single plate, the arrangement of standards, negative controls, and samples was randomized. In keeping with manufacturer instructions, technical replicates remained grouped together within each plate. OD450 was measured using an Epoch Microplate Spectrophotometer (Biotek Instruments, Winooski, VT, United States).
Variation in 5-mC content was modeled as a function of pCO 2 level, exposure time, and their interaction by fitting a linear regression to the data using the lm() function of the R 'stats' package v3.5.1. A full report describing the linear model used to assess changes in 5-mC content is available in Supplementary Table S1. Reported F-statistics from the fitted model were generated using the anova() function of the R 'stats' package. Cliff 's delta estimates, mean effect size estimates, and post hoc Mann-Whitney U-tests comparing 5-mC content between groups were performed with the 'dabest' package for estimation statistics (Ho et al., 2019).
The presence of outliers within the data was assessed by measuring the non-normality, residual, and leverage of each individual replicate fitted to the linear model using Q-Q plots, residual vs. leverage plots, and Cook's distance estimates executed with the R 'stats' base package. Any sample exhibiting nonnormality, high leverage, and high residual was excluded from the linear model, resulting in the removal of a single data point.

RNA Sequencing and Transcriptomic Analyses
RNA was extracted and sequenced from three pooled L. helicina antarctica replicates (10 individuals/pool) per treatment for each of the two sampling time points (n = 18 libraries; 2 µg library −1 ) as previously described (Johnson and Hofmann, 2017). Sequencing was performed at the UC Davis Genome Center on an Illumina HiSeq4000 sequencer using paired-end 100 bp reads. Raw sequencing reads were trimmed using Trimmomatic to remove adapter sequences, low quality base pairs (PHRED < 20), and sequences smaller then 75 bp (Bolger et al., 2014). The quality of trimmed reads was assessed using FastQC v0.11.8 (Andrews, 2018). Trimmed reads were mapped to the L. helicina antarctica reference transcriptome (Johnson and Hofmann, 2016) and counted using RSEM v1.3.1. RSEM was chosen for its accurate quantification of read counts mapped to de novo reference assemblies (Li and Dewey, 2011). Read processing, alignment, and counting was performed with support from the Indiana University Carbonate computing cluster (Stewart et al., 2017). A differential expression analysis was performed using edgeR v3.24.3 (Robinson et al., 2010) set to 'robust' calculation of dispersion and 'robust' fitting of general linear models using the glmQLFit() function of edgeR. Differentially expressed genes (DEGs) were identified using an FDR cutoff < 0.05 and an absolute log 2 foldchange cutoff of > 1.5 . A principal coordinates analysis was also performed in order to cluster RNAseq samples by inputting log 2 -adjusted counts per million (CPM) read counts to the pcoa() function of the R package 'ape' (Paradis et al., 2004). R scripts used to analyze differential gene expression are included in Supplementary Data S3. Enriched gene ontologies (GO) were identified among up-and down-regulated genes using Mann-Whitney U-tests input with signed, −log p-values using the 'Rank Based Gene Ontology Analysis with Adaptive Clustering' method: https:// github.com/z0on/GO_MWU (Wright et al., 2015).
Germline gene body methylation levels were estimated within exons by calculating observed over expected CpG frequencies (CpGOE) across each transcript within the L. helicina antarctica reference transcriptome (Johnson and Hofmann, 2016) using python scripts written by Dimond and colleagues https:// github.com/jldimond/Coral-CpG (Dimond and Roberts, 2016). Means and distributions of transcript CpGOE values were then compared between up-and down-regulated DEGs within treatment groups using permutation tests and Kolmogorov-Smirnov (KS) tests, respectively. A linear model was used to assess variation in transcript CpGOE as a function of foldchange direction, duration of exposure to high pCO 2 , and their interaction using the R package 'lmPerm' v2.1.0, which generates p-values for linear models using permutation tests (Wheeler and Torchiano, 2016). Residuals within this linear model were bimodal, necessitating a statistical approach such as a permutation test that does not assume a normal distribution. Two-sided KS tests were performed using the ks.test() function of the R package 'dgof ' v1.2 (Arnold and Emerson, 2013). The influence of categorical CpGOE bins (bin size = 0.5) on transcriptional variation represented as the CV of logCPM was assessed using the lm() function of R 'stats.'

RESULTS
Genomic Methylation in L. helicina antarctica Under Current and Future pCO 2 Extremes ELISA-based quantification of 5-methylcytosine in L. helicina antarctica revealed significant effects of pCO 2 (F 1,40 = 10.24; p = 0.0027), and the interaction of pCO 2 and time (F 1,40 = 7.74; p = 0.0082) on DNA methylation. Specifically, acute exposure to future pCO 2 extremes induced hypomethylation of the L. helicina antarctica genome followed by a return to normal levels over Frontiers in Marine Science | www.frontiersin.org time. Exposure to current pCO 2 extremes did not induce detectable changes in genomic methylation relative to ambient pCO 2 (Figure 2).
One replicate exhibited 2.14-fold greater 5-mC content than the overall mean. The non-normality, residual, and leverage of this replicate were examined in order to determine its status as an outlier using (i) Cook's distance, (ii) quantile-quantile plots to assess deviation from normality and (iii) residual-leverage plots. The replicate exhibited a substantially high residual and leverage and was thus removed from all analyses (Supplementary Figure S2). A full report describing the linear model used to assess changes in 5-mC content (Supplementary Table S1) and CV values per replicate (Supplementary Figure S3) are available in Supplementary Data S1.

Transcriptomic Responses to Current and Future pCO 2 Extremes
As genomic methylation influences both magnitude and variation in gene expression in multiple invertebrate lineages, we sought to pair genome-wide methylation data with transcriptomes sequenced for early and late timepoints under low, medium, and high pCO 2 . Overall, we found that gene expression in L. helicina antarctica varied as a function of pCO 2 and FIGURE 2 | Differential methylation of the Limacina helicina antarctica genome during acute exposure to future pCO 2 extremes. Mean values (diamonds) and individual replicates (circles) of percent 5-methylcytosine quantified using ELISA-based methods are plotted for 255 µatm (blue), 530 µatm (yellow), and 918 µatm (red) pCO 2 treatments level across time. ± SE is represented by wide, black error bars. ± SD is represented by narrow, colored error bars. exposure time. A principal coordinates analysis of CPM read counts demonstrated that a coordinate axis associated with pCO 2 level and time explained 32.53% of transcriptional variation. A second axis associated with time alone explained 8.4% (Figure 3A).
Limacina helicina antarctica exhibited substantially more differential expression under 918 µatm pCO 2 than under 530 µatm. After 0.5 day of exposure, pteropods in the 530 µatm treatment exhibited only one significant DEG relative to 255 µatm (FDR < 0.05; absolute log 2 FC > 1.5). After 7 days at 530 µatm, 3 transcripts were differentially expressed. Additionally, each DEG in the 530 µatm treatment group was upregulated. By contrast, L. helicina antarctica differentially expressed 6,649 transcripts after 0.5 day of exposure to 918 µatm and 6,815 transcripts after 7 days. 69.91% of DEGs detected at 0.5 day under Frontiers in Marine Science | www.frontiersin.org 918 µatm pCO 2 were downregulated. Downregulated transcripts composed 61.20% of DEGs after 7 days of exposure to 918 µatm ( Figure 3B). 3,607 transcripts were differentially expressed at both 0.5 and 7 days under 918 µatm pCO 2, leaving 3,042 DEGs unique to 0.5 day of exposure and 3,773 DEGs unique to 7 days.
Transcripts that were up-or down-regulated at 0.5 day under 918 µatm pCO 2 were collectively enriched for 88 GO terms while 111 GO terms were enriched among up-and down-regulated transcripts from 7 days of 918 µatm exposure (FDR < 0.01). 81 of these GO terms were commonly enriched among both time points. 56.81% of GO terms enriched at 0.5 day were associated with downregulated transcripts while 54.95% of GO terms enriched at 7 days were associated with downregulation. Ontologies uniquely enriched at 0.5 day pertained to glucose metabolism, protein metabolism, and the oxidative stress response. Ontologies unique to the 7 days timepoint were diverse but included multiple GO terms pertaining to singular functions such as fatty-acid metabolism, antioxidant activity, transcription, trans-membrane transport, and transferase activity for glycosyl and hexosyl groups. Enriched ontologies that were shared between 0.5 and 7 days included functions relating to protein degradation, protein synthesis, cytoskeletal structure, ATP synthesis, and methyltransferase activity (Supplementary Figures S4, S5).
Since methyltransferases are of broad interest in an epigenetic context, we queried enriched GO terms for these subsets of genes. We found that 182 and 174 transcripts associated with methyltransferase activity were significantly enriched among down or upregulated genes at 0.5 and 7 days, respectively, under 918 µatm pCO 2 . On average, enriched methyltransferases were upregulated at both time points. Subontologies of "methyltransferase activity" (GO:0008168) that were enriched among DEGs at both time points included RNA methyltransferases, tRNA methyltransferases, and SAM-dependent methyltransferases: a class that includes enzymes targeting DNA and/or RNA (Supplementary  Figures S4, S5). Histone methyltransferases were not enriched among down or upregulated genes, but some were differentially expressed under future pCO 2 extremes. 1 EZH2 histone-lysine N-methyltransferase was upregulated after 0.5 day exposure while four histone-arginine methyltransferases, 1 UTY histone demethylase, and 4 histone-lysine N-methyltransferases showed downregulation. 1 SETD8 histone-lysine N-methyltransferase was downregulated at both 0.5 and 7 days of exposure to 918 µatm (Supplementary Figures S6, S7).
Full lists of transcripts differentially expressed at 918 µatm pCO 2 relative to 255 µatm are available in Supplementary Data S4, S5. GO terms significantly enriched among sets of up-and down-regulated genes are visualized in Supplementary Figures S4, S5.

Signatures of Gene Body Methylation Among Differentially Expressed Genes
Transcripts that were differentially expressed in response to 918 µatm pCO 2 varied in mean CpGOE ratio depending on fold change direction (F 1,13459 = 85.05; p < 2e-16). CpGOE is inversely related to germline methylation and, when it is applied to transcripts, demonstrates a signature of gene body methylation. Specifically, downregulated genes were enriched with low CpGOE genes at both 0.5 and 7 days ( Figure 4A). A significant interaction between foldchange direction and time was also observed (F 1,13459 = 71.45; p < 2e-16) in which the mean difference in CpGOE between up-and down-regulated DEGs was greater at 0.5 day (Cliff 's delta = 0.1409 ± 0.0330) than 7 days (Cliff 's delta = 0.0292 ± 0.0296).
As ( down-regulated DEGs at 0.5 day (D = 0.19075; p < 1.0e-15) and 7 days (D = 0.15198; p < 1.0e-15) under 918 µatm pCO 2. Density beneath the low CpGOE mode was 2.44-fold greater among downregulated DEGs compared to upregulated DEGs at 0.5 day and 1.84-fold greater after 7 days ( Figure 4A). Differences in the modality of CpGOE among DEGs that were downregulated either at 0.5 and 7 days relative to the whole L. helicina antarctica transcriptome were significant but less pronounced: KS test D values equaled 0.09318-0.06263 for 0.5 and 7 days, respectively (p < 1.0e-15). Cumulative distributions visualizing the modalities of CpGOE in each of the groups are available in Supplementary Figure S8. In contrast to the relationship between genes' CpGOE and differential expression in L. helicina antarctica, a negligible association between CpGOE and transcriptional variation was observed. A linear model reported a significant effect of categorical, binned CpGOE (bin size = 0.5) on the CV of expression of transcripts in which the low CpGOE mode exhibited a higher mean CV (Figure 4B), albeit with a large degree of variation surrounding the mean.

DISCUSSION
Our analyses of DNA methylation and gene expression in L. helicina antarctica showed dynamic changes in both processes following experimental conditioning under future pCO 2 levels predicted for Antarctic coastal waters. Juvenile L. helicina antarctica appeared to mount a strong epigenetic response to future pCO 2 extremes, reducing DNA methylation by 29-56% after brief exposure to high pCO 2 before steadily increasing DNA methylation levels over time. In addition, differential gene expression showed an association with DNA methylation as evidenced by an enrichment of low CpGOE transcripts among genes that were significantly downregulated in response to future pCO 2 extremes. Changes in global DNA methylation and differential expression were absent or negligible following exposure to present-day pCO 2 extremes. Taken together, these observations suggest that DNA methylation plays a role in regulating the cellular responses of L. helicina antarctica to future OA.
Below, we discuss the functional significance of DNA methylation in the context of marine molecular ecology, compare our findings to studies of DNA methylation in other Antarctic organisms, and discuss the available data on epigenetic responses by other species to OA. Lastly, we point to next steps for the study of epigenetics in polar systems, including future work on L. helicina antarctica and other calcifying marine invertebrates.
Hypomethylation Following Acute Exposure to High pCO 2 While L. helicina antarctica has exhibited a poor capacity to acclimate to experimental OA as evidenced by a lack of metabolic compensation and upregulation of stress response pathways following 21 days of acclimation to future pCO 2 levels (Hoshijima et al., 2017;Johnson and Hofmann, 2017), in this study, it appears that juvenile L. helicina antarctica were capable of mounting a significant epigenetic response to high pCO 2 . Our preliminary analyses in L. helicina antarctica have identified dynamic change in levels of DNA methylation in response to high pCO 2 . This outcome might not be expected if we were to stereotype this species as an Antarctic, cold-adapted invertebrate with diminished cellular responses to environmental change. In fact, studies of environmental epigenetics in two Antarctic marine invertebrates, a benthic polychete and an intertidal gastropod, have both documented differential methylation in response to variation in temperature (Marsh and Pasqualone, 2014;Clark et al., 2018).
It is worth noting that one limit to our approach toward quantifying global patterns of DNA methylation is the nonspecific nature of the acquired data. While global levels of methylation provide insight into changes across the whole genome, these measurements are limited in that they (i) provide no indication of functional consequences and (ii) are less sensitive than NGS methods such as bisulfite sequencing. For example, if X methylated CpGs become demethylated in response to a given treatment and are proceeded by methylation at X CpGs that were previously unmethylated, global 5-mC quantification would not detect a change in methylation. For this reason, our results do not rule out the possibility that L. helicina antarctica may differentially methylate its genome in response to presentday pCO 2 extremes. Additionally, the large effect strength of 918 µatm pCO 2 on genomic methylation in L. helicina antarctica is not unprecedented in the context of acute exposure to stress (Rodrigues et al., 2015;Li et al., 2016;Robinson et al., 2019), but the quantity of this effect should be further evaluated with replicate experiments or similar studies.
Although we did detect changes in total methylation and gene expression, our results are presently not sufficient to directly link differential methylation with the expression of specific genes. However, other studies on invertebrates indicate a relationship between gene expression and genomic methylation in response to variation in multiple abiotic factors. Many invertebrate phyla including molluscs predominantly exhibit DNA methylation within gene coding regions (Zemach et al., 2010;Wang et al., 2014;Jeong et al., 2018;Liew et al., 2018). Hypomethylation at intragenic regions is generally associated with reduced gene expression and increased transcriptional variability and/or spuriousness in both vertebrates (Kobayashi et al., 2012;Neri et al., 2017) and invertebrates (Zemach et al., 2010;Dixon et al., 2014Dixon et al., , 2018Gavery and Roberts, 2014;. Many studies have also noted associations between hypomethylation and environmentally-induced or diseaseinduced pathologies spanning a wide range of taxa including plants, vertebrates, and invertebrates for which hypomethylation was demonstrated in both promoters and/or gene bodies (Aina et al., 2004;Pavet et al., 2006;Pogribny et al., 2006;Rusiecki et al., 2008;Luttmer et al., 2013;Xiu et al., 2019). While mechanistic links between hypomethylation and pathology are unclear, there is strong evidence in support of its role in cellular responses to infection and stress.
Despite extensive research on organismal responses to OA, the cellular mechanisms that contribute to pathological vs. adaptive responses to pCO 2 stress remain poorly understood (Melzner et al., 2020). It is plausible that hypomethylation of the L. helicina antarctica genome in response to high pCO 2 is not adaptive. Rather, it may be related to pathological processes associated with chronic stress seen in L. helicina antarctica under OA. The 29-56% reduction in mean 5-mC by L. helicina antarctica in response to pCO 2 stress contrasts observations of differential methylation under OA that have been documented in more pCO 2 -tolerant calcifying invertebrates. Specifically, scleractinian corals and purple sea urchins of differing life history stages respectively exhibited global hypermethylation or a combination of hypo-and hypermethylation at CpGs in response to elevated levels of pCO 2 (Putnam et al., 2016;Liew et al., 2018;Strader et al., 2019;Strader et al., in review). It is unclear to what extent such variation between these taxa is due to phylogenetic differences, life history or to pCO 2 sensitivity, but they remain noteworthy. Further comparative work exploring differences in the variation of DNA methylation under OA in thecosome pteropods relative to more pCO 2 -tolerant calcifying invertebrates could help elucidate mutually exclusive aspects of genomic methylation that contribute to acclimatization and pathology.

Evidence of a Role for DNA Methylation in Differential Expression Under High pCO 2
Limacina helicina antarctica mounted a large transcriptomic response to future pCO 2 extremes that predominantly consisted of downregulated DEGs. This response corroborates our past reporting on differential expression in response to OA among L. helicina antarctica collected from the McMurdo Sound (Johnson and Hofmann, 2017). Since (i) gene expression in molluscs and other invertebrates positively correlates with gene body methylation (Zemach et al., 2010;Gavery and Roberts, 2013;Dixon et al., 2014Dixon et al., , 2018Gavery and Roberts, 2014;, and (ii) DNA methylation occurs primarily within gene bodies among molluscs and other invertebrates, it is possible that hypomethylation of L. helicina antarctica DNA observed in response to high pCO 2 contributed to the differential expression of downregulated genes. Indeed, downregulated DEGs were enriched with low CpGOE transcripts relative to upregulated genes. This demonstrated that downregulated DEGs were likely to have a greater degree of methylation at exons within germline cells, and further, suggested a link between hypomethylation and downregulation across the L. helicina antarctica genome.
Despite evidence of DNA methylation's influence on downregulated transcripts, the timing of changes in 5 m-C and differential expression were not synchronized throughout the duration of the experiment. Mean 5-mC levels in L. helicina antarctica exposed to 918 µatm became comparable to 530 µatm and 255 µatm cultures by 6 days of exposure while the quantity of downregulation in 918 µatm cultures remained high and relatively unchanged even by 7 days. Thus, DNA methylation and differential expression did not correlate over time despite evidence of their association revealed by CpGOE values. DNA methylation and differential expression do not necessarily correlate throughout time when a causal relationship is present (Secco et al., 2015;Pacis et al., 2019), and have been documented to be unpaired during invertebrate development under high pCO 2 (Strader et al., in review). The relationship between DNA methylation and gene expression in molluscs has proved tenuous in some cases. For example, the offspring of diuronexposed oysters have exhibited differential expression and DNA methylation relative to offspring from control parents, but only a small number of DEGs were correlated with differentially methylated CpGs (Rondon et al., 2017). Other reports in oysters have documented stronger associations between gene expression and DNA methylation (Olson and Roberts, 2014). The lack of synchronicity in the timing of differential methylation and differential expression in this experiment does not entirely rule out their mutual influence. Rather, it may suggest a more complex relationship that includes other epigenetic factors.
We also observed interesting patterns of differential expression in genes with various functions related to other epigenetic modifications to DNA. For example, the differential expression of histone modifying enzymes by L. helicina antarctica correlated strongly with changes in global 5-mC over time under high pCO 2 : 1 EZH2 histone-lysine N-methyltransferase was upregulated by L. helicina antarctica after 0.5 day exposure to 918 µatm while 4 histone-arginine methyltransferases, 1 UTY histone demethylase, and 4 histonelysine N-methyltransferases were downregulated. Only one of these genes remained differentially expressed by 7 days of exposure to 918 µatm, demonstrating an early signature of downregulation in histone methyltransferases. In at least some invertebrates and mammals, gene bodies that exhibit higher levels of baseline DNA methylation are also enriched with H3K36me3, histones with lysine 36 trimethylation (Nanty et al., 2011;Baubec et al., 2015). H3K36me3 can chaperone de novo methyltransferases to gene bodies (Baubec et al., 2015). Inversely, the absence of H3K36me3 at gene bodies has been documented to be a predictor of hypomethylation (Hahn et al., 2011). Arginine methylation at H4R3me2 has a similar relationship with de novo methyltransferases but ultimately has a more repressive effect on gene expression (Zhao et al., 2009). Therefore, differential methylation of intragenic regions is associated with chromatin accessibility in certain contexts and both of these processes mediate the influence of intragenic regions on differential expression. Overall, our results suggest that differential methylation of DNA and histones may jointly influence persistent transcriptomic responses by L. helicina antarctica to future pCO 2 extremes. Histone modifications (e.g., histone methylation) and subsequent changes in transcription have been shown to regulate development and responses to environmental stress in molluscs (Fellous et al., 2015(Fellous et al., , 2019Gonzalez-Romero et al., 2017) and other calcifying marine invertebrates (Rodriguez-Casariego et al., 2018). Interactions between chromatin structure and DNA methylation should be explored in greater functional detail in this system in order to substantiate these observations and understand their influence on pervasive and persistent downregulation by L. helicina antarctica under OA.
Pervasive downregulation or "dampening" transcriptomic and/or proteomic responses are commonly observed under high pCO 2 across marine molluscs and other calcifying invertebrates (Todgham and Hofmann, 2009;O'Donnell et al., 2010;Dineshram et al., 2012;Kenkel et al., 2017;De Wit et al., 2018;Kriefall et al., 2018), including some species of Limacina pteropods (Koh et al., 2015;Johnson and Hofmann, 2017). Other genera of shelled pteropods have not exhibited dampening transcriptomic responses to OA (Maas et al., 2015;Moya et al., 2016) and strong time-dependent effects on the proportion of upregulated vs. downregulated genes have been documented in Limacina retroversa in response to high pCO 2 (Maas et al., 2018). Thus, there does not appear to be a generalizable transcriptional response to OA within Pteropoda or the Limacina genus. Studying the potential for phylogenetic variation in epigenetic processes among Pteropods may further explain the contribution of mechanisms such as DNA methylation to cellular responses to OA. Similarly, expanding epigenetic studies of cellular responses to OA across a broader diversity of taxa may reveal mechanistic similarities between evolutionarily distant groups that exhibit similar transcriptional responses to high pCO 2 .

CONCLUSION
Perhaps the most important conclusion we can draw from this experiment is that variation in DNA methylation across the L. helicina antarctica genome was responsive to ocean acidification and likely yielded consequences for cellular functions due to (i) the magnitude of variation, and (ii) evidence of its influence on gene expression. From a polar biology perspective, our results demonstrate that DNA methylation is dynamic in an otherwise stenothermic Antarctic macrozooplankton. Many Antarctic ectotherms have classically been considered to be specialized to an invariable environment, subsequently lacking phenotypic and physiological plasticity necessary for acclimatizing to environmental shifts on par with future climate change (Peck et al., 2004(Peck et al., , 2010(Peck et al., , 2014Beers and Jayasundara, 2015). A growing body of work is beginning to demonstrate, however, that physiological plasticity is retained and sufficient for metabolic recovery under a warmer and/or more acidic ocean in at least some Antarctic ectotherms (Seebacher et al., 2005;Peck et al., 2010;Enzor and Place, 2014;Reed and Thatje, 2015;Huth and Place, 2016a,b;Morley et al., 2016;Enzor et al., 2017;Davis et al., 2018;Hawkins et al., 2018). The responsiveness of DNA methylation in L. helicina antarctica to variation in pCO 2 may be linked to this species' environmental experience in situ. Recent oceanographic observations suggest that seasonal variation in pCO 2 in coastal waters of the Southern Ocean is substantial (Gibson and Trull, 1999;McNeil et al., 2010;Kapsenberg et al., 2015;Negrete-García et al., 2019). Thus, it is likely that L. helicina antarctica is presently faced with low pH conditions that are undersaturated for aragonite calcification during vertical migration.
Lastly, recent modeling efforts have indicated that aragonite undersaturation in the Southern Ocean will appear at quite shallow depths (∼400 m) by the end of the century (Negrete-García et al., 2019) creating a habitat compression or "squeeze" event for calcifying macrozooplankton such as L. helicina antarctica. Thus, quantifying the epigenetic, transcriptional, and physiological plasticity of L. helicina antarctica populations is important in assessing their capacity to respond to future changes in ocean carbonate chemistry. With that said, it remains unclear whether or not dynamic changes in DNA methylation are adaptive for L. helicina populations experiencing high pCO 2 . Future research concerning this ecologically critical species that aims to expand on our findings and past OA experiments will benefit from (i) executing long-term cultures of L. helicina antarctica under future pCO 2 extremes for at least 3-6 weeks and (ii) conducting integrated analyses of metabolic rate, calcification, gene expression, and epigenetic profiling. This 3-6 week range has proven to be a metabolic and transcriptional tipping point for other Antarctic ectotherms acclimating to increases in temperature and pCO 2 (Peck et al., 2010;Enzor and Place, 2014;Enzor et al., 2017;Davis et al., 2018). Quantifying DNA methylation at wholegenome or base-by-base levels across long-term acclimations will expand on dynamic epigenetic changes documented in this study and reveal whether or not they associate with or initiate plastic responses to global change. If DNA methylation or other epigenetic mechanisms ultimately do not drive plastic responses to environmental variation in L. helicina antarctica, this species still represents a valuable comparative system for studying the contribution of epigenetic mechanisms to adaptive vs. pathological responses to abiotic stress and future OA.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the Sequence Read Archive (SRA) BioProject PRJNA576909.

ACKNOWLEDGMENTS
The authors would like to thank the late Dr. Umi Hoshijima for his invaluable assistance throughout the experiment, for collecting in situ pH data reported in this study, and for his consultation regarding these data. The authors acknowledge both Maddie Housh for her assistance in extracting DNA as well as the Crary Laboratory staff at McMurdo Station, Antarctica, for their support. Additionally, other members of the U.S. Antarctic Program, and of Lockheed's Antarctic Support Corporation (ASC) supported field work for this study. The authors also acknowledge the Indiana University Pervasive Technology Institute for providing computational resources on the Carbonate cluster that have contributed to the results reported within this manuscript.