Regulation of gene expression is associated with tolerance of the Arctic copepod Calanus glacialis to CO 2‐acidified sea water

Abstract Ocean acidification is the increase in seawater pCO 2 due to the uptake of atmospheric anthropogenic CO 2, with the largest changes predicted to occur in the Arctic seas. For some marine organisms, this change in pCO 2, and associated decrease in pH, represents a climate change‐related stressor. In this study, we investigated the gene expression patterns of nauplii of the Arctic copepod Calanus glacialis cultured at low pH levels. We have previously shown that organismal‐level performance (development, growth, respiration) of C. glacialis nauplii is unaffected by low pH. Here, we investigated the molecular‐level response to lowered pH in order to elucidate the physiological processes involved in this tolerance. Nauplii from wild‐caught C. glacialis were cultured at four pH levels (8.05, 7.9, 7.7, 7.5). At stage N6, mRNA was extracted and sequenced using RNA‐seq. The physiological functionality of the proteins identified was categorized using Gene Ontology and KEGG pathways. We found that the expression of 151 contigs varied significantly with pH on a continuous scale (93% downregulated with decreasing pH). Gene set enrichment analysis revealed that, of the processes downregulated, many were components of the universal cellular stress response, including DNA repair, redox regulation, protein folding, and proteolysis. Sodium:proton antiporters were among the processes significantly upregulated, indicating that these ion pumps were involved in maintaining cellular pH homeostasis. C. glacialis significantly alters its gene expression at low pH, although they maintain normal larval development. Understanding what confers tolerance to some species will support our ability to predict the effects of future ocean acidification on marine organisms.

The ability of an organism to maintain fitness-related traits, like growth, reproduction, and development, while exposed to a stressful environment by adjusting underlying physiological processes has been referred to as phenotypic buffering (Reusch, 2014). Phenotypic buffering is crucial for tolerance of spatially and temporally variable environments, and organisms that are capable of phenotypic buffering may be able to adapt more easily to long-term environmental change (Donelson, Munday, McCormick, & Pitcher, 2012;Miller, Watson, Donelson, McCormick, & Munday, 2012). However, depending on the mechanism, long-term utilization of phenotypic buffering in a stressful environment may also incur an energetic cost. Phenotypic buffering may involve, among other things, altered gene expression, compensatory feeding, and changes in energy allocation (Sunday et al., 2014).
Recent progress in high-throughput genomic technology has made possible the sequencing of entire transcriptomes (RNA-seq), with the potential to measure expression in all transcribed genes, not limited to a targeted set of genes in a microarray. This has allowed for a universal, nontargeted exploration of gene expression (or transcriptomic) responses to environmental drivers. Investigating the transcriptomic response of species to climate change-related environmental drivers has the potential to provide greater understanding of why and how some species are more tolerant to change than others (Franssen, Gu, & Bergmann, 2011;Reusch, 2014;Seneca & Palumbi, 2015).
Organismal-level measures such as respiration, ingestion, survival, gonad maturation, egg production, and naupliar growth and development are unaffected by low pH (Weydmann, Søreide, Kwasniewski, & Widdicombe, 2012;Hildebrandt, Niehoff, & Sartoris, 2014;Hildebrandt, Sartoris, Schulz, Riebesell, & Niehoff, 2015;Bailey, Thor, & Browman, 2016). In this study, we sought to characterize the molecularlevel response of C. glacialis nauplii to lowered pH to elucidate the physiological processes involved in this tolerance. Understanding the mechanisms of phenotypic buffering that allow some species to tolerate a wide range of pH is important for the field of ocean acidification research. Not only can it help detect low-level stress potentially unnoticed in organismal-level measurements, but it also contributes to understanding how some species, but not others, can tolerate low pH. We quantified whole transcriptome gene expression in the last naupliar stage (N6) of C. glacialis which had been exposed to lowered pH as they developed from eggs (over a month). We hypothesized that gene expression would respond to the pH treatment despite the lack of change in developmental rate, growth, and respiration observed during the same experiment (Bailey et al., 2016) and that these modulations would be related to the maintenance of growth and development.
Furthermore, we expected that the physiological functions affected by pH, as indicated by altered expression of genes coding for functional groups of proteins, would be similar to those affected in other marine organisms at low pH and would, therefore, contribute to an understanding of a generalized molecular response to low pH.

| Collection and exposure
The collection and culturing of the copepods are detailed in Bailey et al. (2016). Briefly, C. glacialis were collected in Rijpfjorden, northeast Svalbard (80°27′54″N, 021°56′63″E) in January 2013 (−1 to 0°C) and held for 50 days prior to the experiment in 40 L flow-through tanks in a cooling container at 2°C. Once females were mature, a pool of 1,900 females was used to inoculate 12 tanks (3 replicate tanks for each of four pCO2 treatments) of 40 L with eggs (7,039 ± 1,178 SE eggs per tank) in turn over 19 days. Females laid eggs in mesh-bottomed buckets suspended in the tanks, which allowed eggs to settle into the tank but retained females. During inoculations, the tanks were filled with ambient sea water; the pH treatment began after the females were removed, with tanks reaching target pH in about 4 hr.
For each treatment, pH was controlled in a 100 L mixing tank by a pH controller (Endress and Hauser, Liquiline CM 442) connected to a pH electrode (Endress and Hauser Orbisint CPS11D glass electrode) calibrated with NIST pH standards (pH NBS ). The controller maintained the target pH by regulating the addition of either CO 2 -enriched sea water (for the low-and mid-pH treatments) or CO 2 -stripped air (for the high-pH treatment) to ambient sea water in the mixing tanks.
The CO 2 -enriched sea water (pH ~5.5) was created by bubbling pure CO 2 continuously into ambient sea water. CO 2 -stripped air was produced by forcing air through a CAS series CO 2 adsorber (Twin Tower Engineering). Sea water from each of the four mixing tanks was pumped into header tanks, each of which fed three 40-L experimental tanks. Inflow of treatment water to the experimental tanks was 10-20 L/min. The pH NBS of each treatment was logged continuously by the pH electrode in each mixing tank. Additionally, water samples (100 mL) from each experimental tank were taken every few days (n = 10 during the experiment) for spectrophotometric measurement of total scale pH (pH T ) with the pH-sensitive indicator dye m-cresol purple (MCP) (Dickson, Sabine, & Christian, 2007). Carbonate chemistry and nutrient concentrations were analyzed collected weekly intervals (n = 8) from one experimental tank of each treatment and preserved with a saturated mercuric chloride solution (Riebesell, Fabry, Hansson, & Gattuso, 2010). Temperature was measured daily in every tank. Methods of spectrophotometric pH T and carbonate chemistry measurements are described in more detail in Bailey et al. (2016). Carbonate chemistry parameters in the experiment are given in Table S9.

| Nauplii sampling
Nauplii were collected from tanks when development had reached stage N6, 35-38 days post-egg laying. The sample was gently sieved and transferred to a petri dish with sea water, where the nauplii were identified to stage by stereoscope. Stage N6 nauplii were transferred to a 10-mL Falcon tube, and RNAlater (Qiagen) was added immediately to preserve the RNA from degradation; the time from collection to preservation was <10 min. Two pools of 10 individuals were collected from each tank. With three replicate tanks per pH, this resulted in six samples from each pH treatment (with five in the pH 7.95 treatment due to sample loss, 23 samples total). The samples were then stored at −80°C for 9 months until extraction.

| RNA-seq
Total RNA was extracted from each of the 23 samples using an and default settings (length fraction = 0.5, similarity fraction = 0.8).
As a proxy of gene expression, the counts of uniquely mapped reads mapping to each contig were compiled for each sample using a custom script from De Wit et al. (2012) with a mapping quality threshold of 20 and a read length threshold of 20 bp.
The proportion of duplicate reads was high. Removal of 110 contigs which matched known C. glacialis ribosomal sequences (1,006 sequences, downloaded from NCBI on 12 April 2016) reduced the number of duplicates, resulting in 45%-93% duplicate reads per sample. As duplicate reads may arise from highly expressed genes and may have biological significance, they were not removed (Parekh, Ziegenhain, Vieth, Enard, & Hellmann, 2015).
To associate physiological function to each of the contigs, the proteins for which the contigs coded were identified by searching the arthropod sequences in NCBI's nonredundant (nr) protein database and the entire UniProt protein database (downloaded on July 6 2016) (e-value cutoff = 1 × 10 −5 ; Blastx 2.3.0+). Hypothetical or predicted protein annotations in nr were excluded by discarding matches with the following words: "hypothetical," "predicted," "unknown," and "putative." The sequences corresponding to known proteins were annotated with biological process and molecular function Gene Ontologies (GO), a classification of physiological function spanning from very broad to specific descriptions for each protein, as well as Kyoto Encylopedia of Genes and Genomes (KEGG) Gene ID using their Uniprot ID and the Uniprot flatfile database (uniprot.org, downloaded on 13 July 2016, http://www. genome.jp/kegg). For each contig with KEGG Gene ID, the KEGG orthologs and their associated KEGG pathways were downloaded from the KEGG REST API using the R package KEGGREST (Tenenbaum, 2016).
KEGG orthologs, and the KEGG pathways they are associated with, are not species specific and are thus designed for cross-species comparisons, allowing their application in nonmodel species like C. glacialis.

| Gene expression analyses
Differentially expressed contigs (DECs), those whose expression varied significantly with pH (treated as a continuous variable), were identified with the R package DESeq2 (version 1.10.1, Love, Huber, & Anders, 2014a). DESeq2 outperforms nearly all other HTS gene expression R packages in detecting significantly differently expressed genes (Seyednasrollah, Laiho, & Elo, 2013;Love, Huber, & Anders, 2014b;Schurch, Schofield, & Gierliński, 2016), taking into account common sources of error such as low sample size, low-count genes and high variability using shrinkage estimators of dispersion and fold change (Love et al., 2014a). Contig expression is discussed here as gene expression.
To determine whether gene expression varied significantly along a continuous pH scale, DESeq2 fits generalized linear models (GLMs) to raw counts of each gene, with a logarithmic link (log2) and negative binomial distribution using a shrunken dispersion estimate (Love et al., 2014a). The GLM fit produces the overall expression strength of the gene and the log fold change (LFC) per unit of pH. LFCs are then shrunk to reduce heteroscedasticity caused by low-count genes showing high LFCs, with low-count, low sample size, or high-dispersion genes shrunk the most. Ward tests with Benjamini-Hochberg procedure are used to adjust for multiple testing in DESeq2. Significant differently expressed contigs were defined as those with |LFC| > 1 and adjusted p-value (false discovery rate, FDR) <0.1.
To visualize the similarity of expression of the DECs across samples, a heatmap was created, with the raw contig counts transformed using the rlog function in DESeq2 and plotted using "aheatmap" in the R package NMF, which centers and scales expression of each contig (0.17.6, Gaujoux & Seoighe, 2010). DESeq2 and heatmap were run on R version 3.2.3 (R Core Team, 2015) and the graphics package ggplot2, version 2.0.0 (Wickham, 2009).
FatiScan is a gene set enrichment analysis (GSEA) tool sensu Khatri, Sirota, and Butte (2012). The first-generation over-representation tools searched for enriched pathways in a subset of genes selected by a threshold of fold change or p-value. However, focusing only on a subset of high-fold-change genes can introduce bias (Dalman, Deeter, Nimishakavi, & Duan, 2012) and distracts from detecting genes that truly affect fitness, which often show mild fold changes (Evans, 2015). In contrast, second-generation GSEA tools use a thresholdfree method that utilizes all genes and their fold changes. They can, therefore, more sensitively detect functional sets of genes (genes in the same GO or pathway) that behave in moderate, but coordinated, manner in response to treatment. Enrichment of GO and KEGG pathways in genes up-and downregulated with pH was analyzed by the FatiScan logistic model using the log fold changes of all genes, with our de novo transcriptome as a reference, and a significance cutoff (α) of FDR-adjusted p-value of .05. GOs with log-odds ratio (LOR) < 0 were taken to be downregulated, and LOR > 0 upregulated.

| RESULTS
The cDNA library sizes differed between treatments and were normalized by DESeq2's coefficient (size factors), which ranged from 0.26 to 3.33, with the size factors of the pH 7.95 treatment being consistently low (Fig. S1). There were 59,353 contigs identified in the de novo transcriptome, with a N50 of 1,019 bp (Supporting Information). BUSCO analysis indicated that the transcriptome was relatively complete: of the 1,066 arthropod orthologs queried, the transcriptome contained 75.4% complete orthologs (71.5% single and 3.9% duplicate copies), with 18.1% fragmented and only 6.5% missing. Of 59,353 contigs, 19,790 (33%) were associated with GO classifications, with an average of 3.54 classifications per contig.
Of these 59,353 contigs, 151 were identified as significantly differently expressed with pH between pH 7.5 and 8.05 (Figure 1). The Numerous cellular processes (1,001 biological process GO terms, 358 molecular function GO terms and 16 KEGG pathways) were significantly downregulated with lowered pH. Clustering these GO terms by their similarity reduced the list to 87 key biological processes, including DNA repair, cell adhesion, growth, cholesterol homeostasis, chromatin modification, development, protein glycosylation, endocytosis, mRNA splicing, regulation of transcription, regulation of protein stability, protein folding, and signaling pathways ( Figure 3; Table S2).
As fewer contigs were significantly upregulated at lowered pH, fewer GO terms (91 biological process GO terms and 39 molecular function GOs) and KEGG pathways (13) were identified as significantly upregulated. The condensed list of biological process GOs that were upregulated at lowered pH include myosin filament organization, locomotion, proteoglycan catabolic process, muscle contraction, transposition, regulation of sodium:proton antiporter, protein deamination, and keratinization ( Figure 4, Table S5). Molecular function GOs upregulated at low pH included motor activity, calmodulin-dependent protein kinase activity, and MAP kinase phosphatase activity (Fig. S4, Table S6). The upregulated KEGG pathways included tight junction, cardiac muscle contraction, pathogenic Escherichia coli infection, ribosome, AMPK signaling pathway, plant-pathogen interaction, and regulation of actin cytoskeleton (Table S7).

| DISCUSSION
In a parallel study on the same experimental cultures, C. glacialis nauplii were able to develop, grow, and respire at normal rates in F I G U R E 1 "Volcano plot" illustrating the selection of genes that were significantly differently expressed with pH in C. glacialis. Plotted are -log-transformed Ward test pvalues for the slope of expression with pH, corrected for multiple tests by the Benjamini-Hochberg method, and presented as false discovery rate (FDR), and -(log2 fold change) of expression over one unit of pH. The log2 fold change (LFC) is negative to indicate that those with a negative value are downregulated as pH is reduced. Significantly differently expressed genes (DECs, in darkest gray) were defined as those with |LFC| > 1 and FDR < 0.1 F I G U R E 2 Heatmap of rlog-transformed gene expression in C. glacialis at different pH treatments, centered, and standardized by gene for each sample (in columns). The color of the cell indicates relative gene expression, with red being higher than the gene's average, and blue lower. Only the 151 significantly differentially expressed genes are shown. The pH treatment and biological replicate (tank) are indicated in shaded red-yellow and gray, respectively, above the heatmap. On the left vertical axis, genes are arranged according to their position in the hierarchical clustering (similarity measure was Pearson correlation coefficient; clustering by complete linkage) CO 2 -acidified sea water down to pH 7.5 (Bailey et al., 2016). On the organismal level, therefore, C. glacialis nauplii were tolerant of pH levels relevant for future ocean acidification scenarios. However, in the present study, significant changes occurred at the transcriptomic level, despite the lack of an observable response on organismallevel vital rates (Bailey et al., 2016). Altered gene expression at low pH indicated regulation of ion exchange activity, DNA transcription, molecular chaperones, and processes involved in locomotor activity in N6 nauplii. This may have allowed C. glacialis to phenotypically buffer the change in pH, thus maintaining normal nauplii development. By measuring the response to low pH at multiple levels of biological organization, both organismal and molecular, we detected physiological changes in response to low pH that were undetected in organismal measures. Additionally, there continues to be considerable F I G U R E 3 Cytoscape network of the 350 most significant downregulated biological process Gene Ontology (GO) terms in C. glacialis, linked by SimRel semantic similarity of the terms by ReViGO. The color of each GO term node (from white to dark red) indicates the adjusted p-value, with more significant terms in darker red. Likewise, the larger the label name, the more significant the adjusted p-value. The size of each node indicates the size of the GO term in the entire UniProt database inconsistency in the measured responses of marine organisms to low pH, even within taxa (Kroeker et al., 2010;Kroeker, Kordas, & Crim, 2013). Therefore, these results on the molecular response of an apparently tolerant species will improve our understanding of the physiological basis of low pH tolerance that is found in some, but not all species.

| Percentage of transcriptome altered
Investigating the expression of single genes identified 151 genes as significantly differently expressed with pH between pH 7.5 and 8.05, representing ~0.25% of the genes in the transcriptome. The de novo transcriptome created in this study was of similar size to that of another C. glacialis assembly (55,562 contigs, Ramos et al., 2015) and . The number of genes that have a significant impact on phenotype is very low, and the percent that modify sensitivity to abiotic stressors is even smaller (0.6%-2% in Brewer's yeast (Saccharomyces cerevisae) deletion strains) (Feder & Walser, 2005;Auesukaree et al., 2009;Evans, 2015). This indicates that even changes in a small percent of the transcriptome can have considerable physiological impact. Acclimation to altered pH over the exposure period may have reduced the number of differentially expressed genes over time, as differential gene expression is often highest immediately after the introduction of the stressor. Huth and Place (2016) found that the degree of altered gene expression in fish exposed to pH and temperature stress decreased with time, being highest 7 days after the start of exposure and decreasing by 80%-90% at 28 and 56 days postexposure start (from 1.0% to 0.1% to 0.2%, respectively). Evans, Chan, Menge, and Hofmann (2013) found a strong reduction in the number of genes differentially expressed with time in purple sea urchin larvae (Strongylocentrotus purpuratus), from 153 at 44 hr of age to only five genes at 96 hr. The acclimation of gene expression can occur over even shorter timescales. Seneca and Palumbi (2015) found many coral genes that were differentially expressed 5 hr after a simulated heat wave but not after 20 hr. The copepods in this study were exposed F I G U R E 4 Cytoscape network of the most significant upregulated biological process Gene Ontology (GO) terms in C. glacialis, linked by SimRel semantic similarity of the terms by ReViGO. Plot components as in Figure 3 for 35-38 days before sampling and therefore may have had larger differences in gene expression between high and low pH immediately after the pH treatments began. Investigating gene expression in N6 nauplii that have been chronically exposed to lowered pH throughout their lives, from the time of hatching through all naupliar stages, reflects an acclimation to low pH rather than the short-term response.
Short-term gene expression changes are unlikely to have a long-term effect on fitness, which is why we chose to focus on the long-term, acclimation response, which more closely reflects how gene expression may respond to near-term climate change with increased pCO 2 in natural systems.
Further, mortality was not tracked in these experiments and is known to be high in broadcast spawned copepod larvae (Kiørboe & Sabatini, 1994). Our measurements, both on the organismal level, as discussed in Bailey et al. (2016), and on the transcriptomic level may therefore reflect the surviving, low pH-tolerant portion of the population. However, this does not change the importance of these data for understanding the transcriptomic basis of tolerance to low pH.

| Regulation of gene set expression
Investigating the coordinated expression of sets of genes revealed that a wide range of biological processes were affected by pH in C. glacialis N6 nauplii. Many of the biological processes downregulated by pH are key components of the evolutionarily conserved, universal cellular stress response (CSR; Kültz, 2003Kültz, , 2005. Cells activate the CSR in response to damage and deformation of macromolecules (proteins, DNA, RNA, and membrane lipids) that occurs as a result of a broad range of stressors (Kültz, 2003). The universal stress response includes the following functions: redox regulation, DNA damage sensing and repair, molecular chaperones, protein degradation, fatty acid/ lipid metabolism, and energy metabolism (Kültz, 2005). In a metaanalysis of 14 studies of transcriptomic responses of oysters (genera  (Table 1).

| Downregulation of stress genes
Interestingly, these CSR-related GO terms and KEGG pathways were downregulated in response to low pH, rather than upregulated. In a short-term response to stress, an upregulation in CSR proteins (or a transcriptomic upregulation of genes coding for them) is expected (Tomanek & Somero, 2000;Kültz, 2003), as has been reported in numerous marine organisms (Lauritano, Procaccini, & Ianora, 2012; Tomanek, 2014; Huth & Place, 2016). However, several studies show both up-and downregulation of CSR-related genes in response to a range of abiotic stressors (Anderson et al., 2015;Goncalves et al., 2016). Downregulation of stress-related genes has been reported in the Sydney rock oyster (Saccostrea glomerata) in response to lowered pH (Goncalves et al., 2016). Downregulation of stress-related genes when exposed to a stressor may also be a characteristic of tolerant populations. In a study comparing wild oysters to a selectively bred line that had developed tolerance to low pH, Goncalves et al. (2016) found that the line that had developed tolerance to low pH reacted to decreased pH by downregulating genes associated with stress response. In contrast, the sensitive wild population upregulated the same genes. Raising wild oysters in low pH water for only three generations produced oysters that responded to low pH stress in the same way as the selectively bred line (selected for fast growth and disease immunity over seven generations), namely downregulation of stress genes (Parker et al., 2012). Similarly, gene expression of an Antarctic fish, whose adaptation to cold waters involves constitutively elevated stress-related genes, shows a downregulation of those genes when exposed to temperature stress (Huth & Place, 2016).
Heat-tolerant strains of Daphnia pulex showed a general downregulation of genes involved in transcription, translation, DNA replication, DNA repair, and core metabolic pathways in response to heat stress, a response not present in heat-sensitive strains (Yampolsky et al., 2014). Downregulation of helicase, involved in DNA replication, was also involved in tolerance to low pH in Pseudocalanus acuspes copepods which had gained some tolerance to low pH after being exposed for three generations (De Wit, Dupont, & Thor, 2015). In this context, the downregulation of the majority of DECs in this study can reasonably be interpreted as an indication that C. glacialis is naturally tolerant of lowered pH, which is consistent with the organismal-level tolerance measured on the same nauplii (Bailey et al., 2016).

| General downregulation
Calanus glacialis nauplii showed an overall downregulation of transcription in the contigs whose expression was significantly related to pH. Of the 151 differentially expressed contigs, 93% were downregulated at low pH. Other studies have also found a general downregulation of proteins and genes in response to low pH stress (Todgham & Hofmann, 2009;Evans & Watson-Wynn, 2014;Dineshram et al., 2016). In a meta-analysis of eight papers evaluating the transcriptomic response of sea urchins to low pH, Evans and Watson-Wynn (2014) found an overwhelming downregulation, with 80% of differentially expressed genes downregulated at low pH (Evans & Watson-Wynn, 2014). Notably, in one paper, 90%-100% of the differentially expressed genes were downregulated (Todgham & Hofmann, 2009). This included cellular stress response genes, metabolism, and apoptosis in urchins exposed to lowered pH, both at pH 7.96 and 7.88. In larvae of the Pacific oyster (Crassostrea gigas), 70% of the differentially expressed proteins were downregulated at low pH, overshadowing the influence of temperature and salinity stressors (Dineshram et al., 2016). They suggested that while the function of general downregulation of protein synthesis and transcription at low pH is unknown, protein synthesis is both energetically costly (Rolfe & Brown, 1997;Wieser & Krumschnabel, 2001;Sokolova, Frederich, Bagwe, Lannig, & Sukhotin, 2012) and produces reactive oxygen species (ROS; Tomanek, 2014). Thus, downregulating protein synthesis at low pH may be a strategy to conserve energy and minimize oxidative damage from ROS. If the downregulation of transcription of some genes observed in C. glacialis nauplii indicates a downregulation of their protein synthesis, this energy conservation strategy may at least partially explain why we did not observe an increase in metabolic rate in the larvae during this experiment (Bailey et al., 2016).
The direction of gene regulation (up-or downregulation) can change with time after onset of stress, with some genes being initially T A B L E 1 Key Gene Ontology (GO) terms and Kyoto Encylopedia of Genes and Genomes (KEGG) pathways regulated significantly with pH in C. glacialis and associated with the universal cellular stress response, as outlined by Kültz (2005) and Anderson et al. (2015). Negative log-odd ratio (LOR) indicates downregulation with lower pH. Number of contigs in the GO term or KEGG pathway (#) and adjusted p-value (p) are also presented Kültz, 2005 Anderson upregulated and then later downregulated in response to heat stress in C. gigas (Meistertzheim, Tanguy, Moraga, Thébault, & Thebault, 2007) or pH and heat stress in the staghorn coral Acropora aspera (Ogawa, Bobeszko, Ainsworth, & Leggat, 2013). In response to stress, Daphnia pulex shows an initial upregulation and then a longer-term downregulation of stress-related genes (Heckmann, Sibly, & Connon, 2008).
This may indicate that the downregulation we measured is a result of the long-term stress exposure prior to sampling. An upregulation of the same stress-related genes may have occurred earlier in the exposure. Depending on the half-life of the proteins they code for, an early upregulation of these stress genes could have produced enough protein that the transcription of more mRNA had been negatively regulated by the time we sampled, thus explaining the general downregulation of many stress-related genes. An exciting future avenue for research would thus be investigation of protein content in combination with transcript-level data over a temporal sampling regime.

| Upregulation
Among the processes upregulated at low pH was the positive regulation of sodium:proton antiporter activity (GO:0032417). The elevated proton concentrations in low pH sea water likely increase the proton concentrations in the extracellular fluid of C. glacialis (Pörtner, 2008;Whiteley, 2011;Melzner, Thomsen, & Koeve, 2013). Aquatic crustaceans maintain pH homeostasis primarily via ion transport (Henry & Wheatly, 1992) and the sodium:proton antiporter is key to maintaining intracellular pH (Casey, Grinstein, & Orlowski, 2010). By removing a proton in exchange for a sodium ion, sodium:proton antiporters directly regulate intracellular pH. Calanus glacialis appears to increase the transcription of genes which upregulate the activity of these antiporters in response to low pH sea water. Other organisms upregulate sodium:proton antiporters in response to low pH water (fish: Catches, Burns, Edwards, & Claiborne, 2006;Huth & Place, 2016) or other ion pumps (sodium:potassium ATPase; Pan, Applebaum, & Manahan, 2015), while corals (Moya et al., 2012;Evans et al., 2013) and urchins (Todgham & Hofmann, 2009;) did not upregulate them at low pH, potentially being able to regulate pH sufficiently with the existing antiporters. Our findings indicate that sodium:proton pumps were upregulated by C. glacialis nauplii in an attempt to maintain cellular pH homeostasis at low pH, although the wide range of altered gene expression at low pH indicates that they did not succeed in fully maintaining homeostasis. Unlike active ion pumps which hydrolyze ATP to move ions across a membrane, sodium:proton antiporters passively exchange ions, a process that does not require the expenditure of metabolic energy (Aharonovitz, Demaurex, Woodside, & Grinstein, 1999). The upregulation of these antiporters, therefore, does not represent an increased energetic cost of pH homeostasis for C. glacialis, indicating that it may be an energetically sustainable component of their response to low pH.
The upregulation of sodium:proton antiporters is also of interest for general copepod physiology. Little is known about the osmoregulation of copepods, and marine copepods are considered osmoconformers (Mauchline, 1998). However, the calanoid copepod Eurytemora affinis has the capability to osmoregulate via the iontransport enzymes V-type H + -ATPase and Na + /K + -ATPase, an ability that allowed it to colonize freshwater from marine environments (Lee, Kiergaard, Gelembiuk, Eads, & Posavi, 2011). Our findings with C. glacialis lend additional support to the osmoregulatory capacity of marine copepods, possibly via sodium:proton antiporters.
Also upregulated were molecular functions related to locomotion, actin-dependent ATPase activity and structural constituents of cuticle and muscle, suggesting altered locomotory muscle use. While we did not measure swimming behavior or escape response speed, larvae of other crustaceans including the American lobster (Homarus americanus) show higher swimming speeds at low pH (Waller, Wahle, Mcveigh, & Fields, 2016). Actinomyosin ATPase is a considerable component of the cellular energy budget (Rolfe & Brown, 1997). If an organism is energy-limited, upregulation of this energy-demanding process may lead to less energy for fitness-related processes such as growth and development, although we did not detect declines in developmental time, body weight, or C:N ratio, or an increase in respiration in the C.

| Transcriptomic expression of tolerance
Understanding the transcriptomic response of an organism to an environmental variable provides information on the extent to which the variable induces a stress response and which physiological functions are induced by that stress response, indicating whether they are components of the universal stress response or specific to that stressor or species. Altered gene expression can indicate stress but also the ability to deal with stressors: Some indicate beneficial responses and some maladaptive.
In evaluating altered gene expression in response to potential stressors, it is important to consider that not all genes have a similar ability to affect fitness. For many genes, the silencing of a single gene produces no phenotypic change at all (Feder & Walser, 2005;Evans, 2015), while other genes have a disproportionate effect on fitness because they are involved in regulatory networks and therefore have the capacity to influence many more proteins than just the one protein that they code for. Additionally, genes which are differentially expressed in response to a stress are not necessarily the same as those that confer tolerance to the same stress (de Nadal, Ammerer, & Posas, 2011;Gibney, Lu, Caudy, Hess, & Botstein, 2013;Giaever & Nislow, 2014;De Wit et al., 2015;Evans, 2015). Two types of genes that usually are not differentially expressed under stress, but whose expression have far-reaching effects that can confer tolerance to that stress, are those involved in chromatin modification and cell signaling (Gibney et al., 2013;Jarolim, Ayer, & Pillay, 2013;Evans, 2015). Interestingly, in C. glacialis, genes controlling both chromatin modification (GO:0016568) and intracellular signaling transduction (GO:0035556) were significantly affected by pH, being downregulated in lower pH.
As both chromatin and signaling pathways have the potential to regulate a large number of other genes and proteins, this indicates that large physiological changes may have been altered in C. glacialis via these transcriptional changes. Chromatin modification affects DNA replication, transcription, and repair and is a critical component of an organism's response to stress (Smith & Workman, 2012) and acclimation to their environment (van Zanten, Tessadori, Peeters, & Fransz, 2013). Changes to chromatin can also be passed on to offspring via epigenetic inheritance, thus potentially affecting the response of future generations to the environment (Campos, Stafford, & Reinberg, 2014). Indeed, epigenetic acclimation to stressors may be possible in copepods. Transgenerational effects (either adaptation or epigenetic acclimation) in the copepod Pseudocalanus acuspes alleviated the negative effect of low pH after two generations in low pH . While the temporal duration of the chromatin modifications in C. glacialis at low pH, or their ability to be epigenetically passed on to the next generation, is unknown, they may play a part in both acclimatization within a single lifetime and over several generations.
The degree of altered gene expression during stress does not necessarily correlate with a reduction in performance or fitness. Activation of the cellular stress response immediately after the onset of the stress can help protect cells against protein unfolding, DNA damage, and oxidative damage (Kültz, 2003(Kültz, , 2005. Differential gene expression can therefore indicate adaptive responses and contribute to tolerance. Smolina et al. (2015) suggested that a strong transcriptomic response to stress underpins tolerance at the organismal level (i.e., respiration, feeding, mortality). They showed that C. finmarchicus strongly altered their gene expression during heat stress, while C. glacialis exhibited a near lack of response, altering the expression of only one gene transcript. As a result, C. finmarchicus showed mild changes in feeding and mortality at increased temperature, while C. glacialis was strongly affected, reducing grazing and showing higher mortality rates than C. finmarchicus under heat stress. In contrast to Smolina et al.'s (2015) findings on C. glacialis' mild gene expression response to increased temperature, we found that C. glacialis exhibited a significant change in gene expression in response to low pH. Similar to their findings, we found that significant regulation of gene expression was associated with tolerance to stress. The gene expression regulation in C. glacialis nauplii is associated with tolerance of low pH for a range of organismallevel vital rates (Bailey et al., 2016).
In evaluating the implication of transcriptomic regulation on tolerance, it is important to consider the time course of changes in gene expression. With the emergence of transcriptomics studies in investigating organismal responses to environmental stressors, the concept of "transcriptomic resilience" has recently been introduced (Franssen et al., 2011). They propose that the quick return to normal gene expression after exposure to a stressor is a sign that an organism can tolerate environmental changes. Seneca and Palumbi (2015) support this theory, showing that gene expression in corals from a tolerant population returned to control levels to a greater extent than in corals from a sensitive population, in response to heat stress. Those individuals whose gene expression returned to control levels exhibited less bleaching than those whose gene expression remained altered.
Tracking the temporal signal of differential gene expression is therefore of utmost importance when comparing the stress tolerance of different species or the effect of different stressors. To achieve this, we suggest that future studies on gene expression responses to stressors collect samples over time, starting immediately after exposure initiation, in order to track the temporal change in gene expression.
The concentration of specific proteins is not only the result of transcription, but also the balance between mRNA degradation, translation, and protein degradation (Schwanhäusser, Busse, & Li, 2011). While stress can alter translation (Picard, Loubière, Girbal, & Cocaign-Bousquet, 2013), mRNA degradation (Huch & Nissan, 2014), protein degradation (Flick & Kaiser, 2012), and protein activity (Krebs & Holbrook, 2001;Pan et al., 2015), transcriptomic regulation nonetheless can show a clear pattern of response to environmental stressors and, thus, indicate the involvement of genes and physiological processes in the organism's response (de Nadal et al., 2011;Evans & Hofmann, 2012;Evans & Watson-Wynn, 2014;Anderson et al., 2015). With six measurements at each pH (five in the pH 7.95 treatment) and samples consisting of pooled individuals, we were able to conservatively select only the consistent signals in gene expression (Schurch et al., 2016). Furthermore, few papers examine differential gene expression over >2 pH treatments or use pH as a continuous variable (i.e., Barshis et al., 2013;Seneca & Palumbi, 2015;Huth & Place, 2016). By treating pH as a continuous variable (as opposed to a factor) and collecting data from four pHs, we narrowed the focus of this study to genes differentially expressed consistently over the investigated pH range, rather than those reacting to specific treatments, and reduced the chance of false positives.
Altered gene expression allowed C. glacialis to phenotypically buffer the change in pH, thereby maintaining normal nauplii development.
The observed downregulation of stress-related genes in response to stress appears to be a characteristic of stress-tolerant populations and species. Additionally, the upregulation of energy-neutral ion exchangers may contribute to C. glacialis' response to low pH without adding additional cost. In predicting the effects of ocean acidification on marine organisms, understanding the molecular mechanisms that confer tolerance to some species will support our ability to predict the effects of future ocean acidification on marine organisms.

ACKNOWLEDGMENTS
The authors would like to acknowledge support from Science for Life Laboratory, the National Genomics Infrastructure, NGI, and Uppmax for providing assistance in massive parallel sequencing and computational infrastructure. We thank NOTUR (UNINETT Sigma2) for com-

CONFLICT OF INTEREST
None declared.

DATA ACCESSIBILITY
All short Illumina reads were submitted to NCBI SRA (accession numbers SAMN05990738-SAMN05990760; BioProject ID: PRJNA352656). The transcriptome assembly is provided in the Supporting Information.