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₂ 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 winter-summer transition in the McMurdo Sound is characterized by low phytoplankton abundance (Foster et al., 1987).

Seawater pCO₂ was manipulated using three reservoir containers held at different pCO₂ levels, each connected to three replicate culturing tanks of the volume and flow rate described above. pCO₂ was controlled inside of reservoir tanks by injecting mixtures of filtered, CO₂-scrubbed, dry air and pure CO₂ with seawater. Pure air and CO₂ were mixed using SmartTrak^® 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₂, pH, and carbonate chemistry were measured as previously described (Johnson and Hofmann, 2017). Low, intermediate, and high pCO₂ 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.

DNA was extracted from pooled replicates containing 10 juvenile L. helicina antarctica sampled from each pCO₂ 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₂ 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₂ 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₂ 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 non-normality, high leverage, and high residual was excluded from the linear model, resulting in the removal of a single data point.

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₂ foldchange cutoff of > 1.5 (Chen et al., 2016). A principal coordinates analysis was also performed in order to cluster RNAseq samples by inputting log₂-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₂, 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.’

ELISA-based quantification of 5-methylcytosine in L. helicina antarctica revealed significant effects of pCO₂ (F₁,₄₀ = 10.24; p = 0.0027), and the interaction of pCO₂ and time (F₁,₄₀ = 7.74; p = 0.0082) on DNA methylation. Specifically, acute exposure to future pCO₂ extremes induced hypomethylation of the L. helicina antarctica genome followed by a return to normal levels over time. Exposure to current pCO₂ extremes did not induce detectable changes in genomic methylation relative to ambient pCO₂ (Figure 2).

Individuals exposed to 918 μatm pCO₂ (future extremes) for 1 day exhibited a 0.42% (0.29–0.56 95% CI) reduction in 5-mC content relative to 255 μatm (Mann–Whitney p = 0.020). For context, the mean 5-mC content of the L. helicina antarctica genome equaled 0.944 ± 0.17% when animals exposed to 918 μatm pCO₂ were excluded from the mean. 5-mC in 918 μatm cultures remained 0.19% lower than the 255 μatm treatment after 3 days of exposure (Mann–Whitney p = 0.060) before converging at a difference of 0.0039% after 6 days. L. helicina antarctica showed no differences in 5-mC content between 530 and 255 μatm pCO₂ treatments at any timepoint (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.

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₂. Overall, we found that gene expression in L. helicina antarctica varied as a function of pCO₂ and exposure time. A principal coordinates analysis of CPM read counts demonstrated that a coordinate axis associated with pCO₂ 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₂ 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₂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 918 μatm pCO₂ 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₂, 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₂ 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₂. On average, enriched methyltransferases were upregulated at both time points. Sub-ontologies 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₂ 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₂ 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.

Transcripts that were differentially expressed in response to 918 μatm pCO₂ varied in mean CpGOE ratio depending on fold change direction (F₁,₁₃₄₅₉ = 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₁,₁₃₄₅₉ = 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 (i) the distribution of CpGOE in the L. helicina antarctica transcriptome is bimodal and (ii) the maximum and minimum percentiles of CpGOE distributions were equivalent between up- and down-regulated transcripts, changes in the densities of CpGOE distributions among DEGs served as a more accurate measurement of the frequency of high and low CpGOE values. A Kolmogorov–Smirnov (KS) test demonstrated significant differences in the modality of CpGOE values between up- and 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.

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₂ 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₂. Our preliminary analyses in L. helicina antarctica have identified dynamic change in levels of DNA methylation in response to high pCO₂. 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 non-specific 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 present-day pCO₂ extremes. Additionally, the large effect strength of 918 μatm pCO₂ 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; Li 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., 2014, 2018; Gavery and Roberts, 2014; Li et al., 2018). Many studies have also noted associations between hypomethylation and environmentally-induced or disease-induced 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₂ stress remain poorly understood (Melzner et al., 2020). It is plausible that hypomethylation of the L. helicina antarctica genome in response to high pCO₂ 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₂ stress contrasts observations of differential methylation under OA that have been documented in more pCO₂-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₂ (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₂ 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₂-tolerant calcifying invertebrates could help elucidate mutually exclusive aspects of genomic methylation that contribute to acclimatization and pathology.

Limacina helicina antarctica mounted a large transcriptomic response to future pCO₂ 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., 2014, 2018; Gavery and Roberts, 2014; Li et al., 2018), 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₂ 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₂ (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 diuron-exposed 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₂: 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 histone-lysine 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₂ 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, 2019; Gonzalez-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₂ 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₂ (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₂.