Gene expression and epigenetic responses of the marine Cladoceran, Evadne nordmanni, and the copepod, Acartia clausi, to elevated CO2

Abstract Characterizing the capacity of marine organisms to adapt to climate change related drivers (e.g., pCO2 and temperature), and the possible rate of this adaptation, is required to assess their resilience (or lack thereof) to these drivers. Several studies have hypothesized that epigenetic markers such as DNA methylation, histone modifications and noncoding RNAs, act as drivers of adaptation in marine organisms, especially corals. However, this hypothesis has not been tested in zooplankton, a keystone organism in marine food webs. The objective of this study is to test the hypothesis that acute ocean acidification (OA) exposure alters DNA methylation in two zooplanktonic species—copepods (Acartia clausii) and cladocerans (Evadne nordmanii). We exposed these two species to near‐future OA conditions (400 and 900 ppm pCO2) for 24 h and assessed transcriptional and DNA methylation patterns using RNA sequencing and Reduced Representation Bisulfite Sequencing (RRBS). OA exposure caused differential expression of genes associated with energy metabolism, cytoskeletal and extracellular matrix functions, hypoxia and one‐carbon metabolism. Similarly, OA exposure also caused altered DNA methylation patterns in both species but the effect of these changes on gene expression and physiological effects remains to be determined. The results from this study form the basis for studies investigating the potential role of epigenetic mechanisms in OA induced phenotypic plasticity and/or adaptive responses in zooplanktonic organisms.

The impacts of ocean acidification (OA) on marine organisms, are complex and species-specific (Fabry et al., 2008;Kroeker et al., 2013;Kurihara, 2008;Wittmann & Pörtner, 2013). As genomic analysis techniques increase in efficiency, transcriptomic analyses have become important metrics for quantifying the expression of stressrelated genes (Evans & Hofmann, 2012), and enable the examination of a broad range of genetic responses to environmental change (Harms et al., 2014;Todgham & Hofmann, 2009). These techniques increase our understanding of the scope of the organismal response and have the potential to detect molecular compensation for environmental stress that may otherwise go undetected using more traditional physiological studies (Gracey, 2007). For example, in the copepod Calanus glacialis exposure to high pCO 2 -low pH significantly alters genes associated with cellular stress response, oxidative stress and ion transporters, suggesting an important role for these genes in pH homeostasis (Bailey et al., 2017). Similar transcriptional responses have been reported in many other species (Strader et al., 2020). OA conditions alter DNA methylation patterns in several species of invertebrates (Downey-Wall et al., 2020;Liew et al., 2018;Putnam et al., 2016). However, there is no experimental evidence demonstrating a direct association between molecular changes and physiological/morphological phenotypes.
The regulation of eukaryotic gene expression occurs through many pathways, including but not limited to epigenetic mechanisms such as DNA regulation and histone modifications (Gibney & Nolan, 2010).
Methylation of genomic DNA provides a mechanism to modulate transcription by influencing the binding of regulatory factors to regulatory elements (Medvedeva et al., 2014). In vertebrates, DNA methylation is involved in the regulation of gene expression during cellular differentiation and development (Law & Jacobsen, 2010). Methylation of CpG dinucleotides, particularly in the promoter regions immediately upstream of the transcription start site, play a major role in the regulation of gene expression (Suzuki & Bird, 2008). More recently, it has become clear that DNA methylation in other genomic regions, such as intergenic and gene body methylation, also play important roles in gene regulation (Jeziorska et al., 2017;Zhou et al., 2020). In insects, DNA methylation is restricted to the transcribed regions and plays an important role in behavioral plasticity and social behavior (Yan et al., 2014). Research on DNA methylation is rare in crustaceans and has focused largely on the freshwater cladoceran, Daphnia spp. (Kvist et al., 2020;Lindeman et al., 2019).
The objective of this study is to investigate the effect of acute exposure to OA on transcriptional and epigenetic (DNA methylation) patterns in two marine planktonic crustaceans, the copepod Acartia clausi and the Cladoceran Evadne nordmanni. While copepods reproduce sexually, Cladocerans have both sexual and parthenogenic reproductive pathways. Cladocerans are widely used as models to study the evolutionary basis of phenotypic plasticity because they reproduce clonally (asexually) and sexually, which offers a unique opportunity to assess the relative contributions of the epigenetic (in clonal populations) and genetic (in sexually reproducing populations) mechanisms underlying adaptation to environmental drivers and their molecular basis (Harris et al., 2012;Lindeman et al., 2019;Toyota et al., 2019). By comparing the responses of copepods and cladocerans (asexually vs. sexually reproducing populations) we hope to establish the relative roles of genetic and epigenetic mechanisms in determining the adaptation capacity of marine populations to pCO 2 . The results from this study provide a basis for investigating the possibility for rapid adaptive responses of planktonic marine organisms to climate change. Until recently, studies addressing the epigenetic basis of adaptation have been conducted mainly in sessile animals such as corals and mollusks. The system used to augment the CO 2 levels of the seawater is described in Runge et al (2019). Briefly, seawater is pumped from the Bjørnafjord at a depth of 160 m up to the laboratory facilities, where it is sand-filtered then passed through a 20 mm Arcal disk filter. This seawater input is then split between a set of 100 L mixing tanks and a stock solution tank, that is, bubbled with CO 2 to maintain a pH of 5.8. Dosing pumps (Iwaki Inc.), controlled by feedback from pH electrodes and controllers (Endress and Hauser, Liquiline CM 442), then add the low-pH stock solution to the mixing tanks to create seawater at treatment pCO 2 levels (Table 1).

| Animal collection and experimental design
Animals were placed in 250 ml culture flasks filled with water from this system and held in a temperature-controlled environmental chamber at 12.5°C (±0.1°C) for 24 h. Experiments were run with four replicate chambers with ~20 E. nordmanni or ~40 A. clausi in each chamber. Full carbonate chemistry of the water was measured prior to exposure and pH was measured after exposure. Given the sealed environments and limited exposure time, electrode measurements at the end of 24 h showed that pH never increased by more than 0.068 units in any replicate, averaging a 0.027 unit increase from the initial values. At the end of the experiment, animals were filtered from each chamber and pipetted into 3 ml cryotubes and frozen at −80°C until analyzed for gene expression and DNA methylation profiles.

| Carbonate chemistry
Carbonate chemistry was calculated from the total alkalinity (AT; measured by titration) and from pH. The pH (total) was measured spectrophotometrically (Hitachi U-2900 dual-beam) using the pHsensitive indicator dye m-cresol purple (Sigma-Aldrich) following SOP (standard operating procedure 6b: (Dickson et al., 2007)). Samples of seawater were collected in 20 ml scintillation vials (leaving no head space) from all experimental vessels and held in a dark, 25°C water bath for temperature equilibration. pH was always measured within 3-5 h of sample collection. To make each pH measurement, 10 ml of each sample was slowly pipetted into two quartz cuvettes with a 5 cm path length (a modification of the 10 cm path length in SOP 6b). The cuvettes were sealed with a Teflon cover, and held at 25°C in the temperature-controlled chamber of the spectrophotometer.
M-cresol purple (10 μl) was added to the sample cuvette, while the second cuvette served as a reference. Absorbance was measured at 578 nm (A1), 434 nm (A2), and 730 nm (background). We used equations in section 8.3 of SOP 6b to correct A1/A2 for the addition of dye. The pK2 and final pH value was determined from Liu et al. (2011, Equation 18). Carbonate chemistry was determined from pH, total alkalinity (AT), temperature, salinity, and nutrients (phosphate and silicate). AT was analyzed by potentiometric titration (Dickson et al., 2007) in an open cell with 0.1 M HCl using a Hach AT1122 [SS1] automatic titrator (Loveland, CO, USA). Certified reference material provided by Andrew Dickson (Scripps Institution of Oceanography, San Diego, USA) was used to calibrate AT measurements. An additional sample (20 ml) was collected and stored in HDPE bottles with HDPE caps with 200 μl of chloroform added. These samples were then analyzed for silica, phosphorus, and nitrogen. Additional carbonate chemistry parameters were calculated using CO2SYS2.3 (Lewis et al., 1998) with the standard set of carbonate system equations and constants of (Mehrbach et al., 1973) after applying the refit of (Dickson & Millero, 1987).

| Transcriptome profiling
Unstranded RNAseq libraries for both the species were prepared using the Illumina TruSeq total RNA library prep kit and 50 bp singleends sequencing on the HiSeq2000 platform were performed at the Tufts University Core Facility. Raw data files were preprocessed by trimming the reads using Trimmomatic (Bolger et al., 2014), removing the low quality reads (Phred score <35) and assessed for quality using FastQC (Andrews, 2010). Preprocessed reads were concatenated and de novo transcriptome assembled using Trinity (v.2.8.6, (Haas et al., 2013)) ( Figure 1a). We used default parameters except for minimum contig length (set to 300) and normalize maximum read coverage (set to 50). Functional annotation was performed by first predicting the coding regions of the transcripts using TransDecoder (v3.0.0) (Haas et al., 2013). The predicted transcripts were annotated using BLASTx (Altschul et al., 1997). The resulting statistically significant BLAST annotations were used for Gene Ontology (GO) classification, a system for hierarchically classifying genes based on their biological process and molecular function (Ashburner et al., 2000). Differential gene expression was carried out by mapping the reads to the de novo assembly using kallisto (Bray et (Robinson et al., 2010). GO analysis of differentially expressed genes (DEGs) was done using gProfiler. Bonferroni correction for multiple testing (p-value < .05) was used while determining the fold enrichment. To understand the relationship between GO terms, Directed Acyclic Graphs of significantly enriched GO terms were drawn using GOView (webgestalt.org/GOView). Raw data files have been deposited in NCBI BioProject (PRJNA91434).

| Quantification of global DNA methylation levels
Global cytosine levels were determined using established high performance liquid chromatography (HPLC) protocols (Ramsahoye, 2002) with some modifications. Briefly, genomic DNA was hydrolyzed with a combination of RNase A and RNase T1 followed by ethanol precipitation. DNA was digested to nucleosides as described previously (Hashimoto et al., 2015) and separated on a 6490 Triple Quad LC-MS with UV absorbance detector (1290 Infinity UV detector, 6490 Triple Quad Mass detector, Agilent, Santa Clara, CA) equipped with an XSelect™ HSS T3 column (2.1 × 100 mm, 2.5 μm, Waters, Milford, MA). Adaptor-ligated fragments of 150-250 bp and 250-350 bp in size were recovered from a 2.5% NuSieve 1:1 agarose gel (Zymoclean™ Gel DNA Recovery Kit). The fragments were then bisulfite-treated using the EZ DNA Methylation-Lightning™ Kit. Preparative-scale PCR was performed and the resulting products were purified and 50 bp paired end (PE) sequencing was performed on an Illumina HiSeq2500 platform. Sequence reads from eRRBS libraries were identified using standard Illumina base-calling software. Raw reads were preprocessed using TrimGalore and aligned to the genome using RefFreeDMA (Klughammer et al., 2015). RefFreeDMA was designed for conducting differential DNA methylation analysis with deduced reference genome, thus allowing DNA methylation profiling in organisms without a reference genome (Figure 1b). We conducted DNA methylation analysis only in the CpG context. We did not conduct the bisulfiteBlast step in RefFreeDMA pipeline to verify species annotations due to several issues in modifying the provided script for our experimental dataset. Raw data files have been deposited in NCBI BioProject (PRJNA780490).

| Bisulfite PCR (BS-PCR)
Methylation analysis of CpG islands was performed by BS-PCR. A 50 μl PCR reaction was carried out in 1X PCR buffer, 5 mM MgCl 2 , 1 mM dNTP mix, 1 unit of Taq polymerase, 50 pmol each of the forward and reverse primers, and ~50 ng of bisulfite-treated genomic DNA. BS-PCR primers were designed using the sense strand of the bisulfite-converted DNA. PCR cycling conditions were 94°C for 10 min, followed by 40 cycles of (94°C for 30 s, 55°C for 30 s and 72°C for 30 s), and a final cycle of 72°C for 8 min. PCR products were electrophoresed on 1% agarose gels, bands excised and gel-extracted using the Gene Clean II kit (MP Biomedical, CA). Purified PCR products were cloned using the pGEM-Teasy cloning kit (Promega, MI) as per the manufacturer's protocol. Mini-preps were prepared using the Pure Yield plasmid miniprep kit (Promega, MI). The primer sequences are provided in Table S1. For each sample, a minimum of 5 clones were sequenced. BS-PCR together with sequencing of several clones provides allele-specific methylation profiles.

| RE SULTS
Calculated pCO 2 values in the low pCO 2 exposures ranged from 455 (14 SD) to 552 (31 SD) μatm (Table 1)  We compared the DEGs observed in response to high pCO 2 exposure in both species and observed 49 DEGs to be commonly expressed. Among them, 16 genes were upregulated and 33 downregulated. GO analysis suggests that many of these genes fall under the GO term one carbon metabolic process (GO:0006730), oxidative phosphorylation (GO:0006119), and NADPH regeneration (GO:0006740). We observed five DEGs related to the GO term-one carbon metabolic process in both species (Figure 3).   conditions are highly variable between these studies, there are some consistent transcriptional responses observed in multiple species.

| D ISCUSS I ON
One of the well documented changes includes differential expression of genes associated with metabolism. We observed similar responses in both species. This is not surprising given the fact that elevated carbon dioxide causes metabolic depression (Michaelidis et al., 2005;Pörtner et al., 1998;Reipschläger & Pörtner, 1996). This mainly results from a decrease in extracellular pH and compensatory mechanisms exist to re-establish the acid-base regulation. We observed differential expression of genes encoding ion transporters, suggesting that mechanisms associated with metabolic depression are highly conserved.
Another widely demonstrated effect of high CO 2 exposure is al-   (Donelson et al., 2012;Sunday et al., 2014). However, these results are species-specific. It is widely established that persistent changes in physiology or gene expression, that is, induced by a previous exposure to stressors are considered an "epigenetic" memory . This memory is based on DNA methylation and histone modifications that alter chroma- the GO pathways that was enriched in the RNAseq dataset is the one-carbon metabolism pathway. The substrates from one-carbon metabolism play an important role in the maintenance of cellular nutritional status by converting nutrients (e.g., glucose, amino acids) into metabolites that feed into diverse biological functions, including cellular biosynthesis, maintaining cellular redox status. In addition, it provides substrates involved in the regulation of protein and nucleic acid methylation (Locasale, 2013). The observed CO 2 exposure related differences in one-carbon metabolic pathway suggest that cellular metabolism is impacted. This could potentially affect epigenetic regulation of gene expression, nucleic acid biosynthesis and metabolic disturbances, something that should be investigated in future studies.
Global 5-methyl cytosine profiling of DNA using two different methods did not reveal any significant changes in response to high pCO 2 exposure in either of the species. However, genome wide profiling identified several differentially methylated regions with two-thirds of them hypomethylated in response to OA. While these results suggests that methylation or demethylation may have an important role in the organisms response to OA, the lack of a well annotated genomic resource for these species makes assigning differentially methylated regions to specific genes impossible. Future work in this arena is certainly warranted.
In recent years, the belief that the genetic code is the sole basis for biological inheritance has been challenged by the discovery of trans-generational epigenetic inheritance. Through epigenetics, environmentally induced phenotypes can persist for several generations, due to the transmission of molecular factors that determine how DNA is read and expressed (Bonduriansky & Day, 2009;Jablonka & Raz, 2009;Verhoeven et al., 2016). While epigenetic inheritance is well documented (Verhoeven et al., 2016), the adaptive significance, if any, of such a complementary inheritance system remains enigmatic. Since it constitutes the inheritance of an environmentally induced phenotype, its adaptive value should depend upon whether environments are predictable across generations. The use of clonal species such as cladocerans provides a promising avenue for determining the degree of heritability and the longevity of epigenetic changes in subsequent generations.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
High throughput sequencing raw data (RNA sequencing and DNA methylation) has been deposited in NCBI BioProject (Accession number: PRJNA780490).