Acclimatory gene expression of primed clams enhances robustness to elevated pCO2

Sublethal exposure to environmental challenges may enhance ability to cope with chronic or repeated change, a process known as priming. In a previous study, pre‐exposure to seawater enriched with pCO2 improved growth and reduced antioxidant capacity of juvenile Pacific geoduck Panopea generosa clams, suggesting that transcriptional shifts may drive phenotypic modifications post‐priming. To this end, juvenile clams were sampled and TagSeq gene expression data were analysed after (i) a 110‐day acclimation under ambient (921 μatm, naïve) and moderately elevated pCO2 (2870 μatm, pre‐exposed); then following (ii) a second 7‐day exposure to three pCO2 treatments (ambient: 754 μatm; moderately elevated: 2750 μatm; severely elevated: 4940 μatm), a 7‐day return to ambient pCO2 and a third 7‐day exposure to two pCO2 treatments (ambient: 967 μatm; moderately elevated: 3030 μatm). Pre‐exposed geoducks frontloaded genes for stress and apoptosis/innate immune response, homeostatic processes, protein degradation and transcriptional modifiers. Pre‐exposed geoducks were also responsive to subsequent encounters, with gene sets enriched for mitochondrial recycling and immune defence under elevated pCO2 and energy metabolism and biosynthesis under ambient recovery. In contrast, gene sets with higher expression in naïve clams were enriched for fatty‐acid degradation and glutathione components, suggesting naïve clams could be depleting endogenous fuels, with unsustainable energetic requirements if changes in carbonate chemistry persist. Collectively, our transcriptomic data indicate that pCO2 priming during post‐larval periods could, via gene expression regulation, enhance robustness in bivalves to environmental change. Such priming approaches may be beneficial for aquaculture, as seafood demand intensifies concurrent with increasing climate change in marine systems.

Organismal environmental resistance depends on integration of predictable environmental cues into acclimatory phenotypes (Ghalambor et al., 2007;Snell-Rood et al., 2018). Although larvae of marine metazoans are highly susceptible to changes in the surrounding environment, early life presents an ideal window for developmental acclimatization due to the importance of environmental information in setting the stage for subsequent phenotypic outcomes (Burton & Metcalfe, 2014;Fawcett & Frankenhuis, 2015). Environmental variation (both spatial and temporal) shapes phenotypes (Dowd et al., 2015) and numerous studies support an acclimatory capacity for marine invertebrates to cope with intermittent periods of thermal stress (Hraoui et al., 2021) and elevated pCO 2 within a generation (Détrée & Gallardo-Escárate, 2018;Gurr et al., 2020;Li et al., 2020;Suckling et al., 2015), as well as across a generation (Goncalves et al., 2016;Parker et al., 2015). Thus, the timing and magnitude of environmental change is likely to have a joint effect on plasticity . Carryover effects of environmental history have physiological (Espinel-Velasco et al., 2021;Parker et al., 2012), ecological (Costantini, 2014;Hettinger et al., 2013) and evolutionary implications (Thomsen et al., 2017). In light of this, it is essential to understand how intermittent or repeated environmental signals, such as the challenges posed from climate change, are transduced to elicit acclimatory mechanisms.
Gene expression regulation is key to homeostasis and phenotypic plasticity. Both constitutive and inducible gene expression responses can enhance acclimatory capacity (Barshis et al., 2013;Georgoulis et al., 2021). Transcriptome profiling of clams and oysters exposed to changing environmental conditions has identified differential regulation of mitochondrial complexes, antioxidants and proteins related to lipid degradation (Chapman et al., 2011;Goncalves et al., 2017;López-Landavery et al., 2021;Teng et al., 2021), indicating that external abiotic conditions can affect metabolism and shift substrates for bioenergetics. Furthermore, enhanced constitutive expression, or gene frontloading, is a proposed mechanism to cope with challenging but predictable environmental signals (Barshis et al., 2013). For instance, limpets (Lottia sp.) occupying the high intertidal zone transcribe heat-shock proteins at a higher level relative to low-intertidal individuals, suggesting an acclimated response (Dong et al., 2008). Thus, transcriptomics provides a broad and sensitive means to assess the role of gene expression in environmental priming and the molecular underpinnings of important economic traits in aquaculture species (Chandhini & Kumar, 2019).
Transcriptome profiles of geoduck clams reared under OA conditions found regulation of energy production and acid/base homeostasis in umbonate-stage larvae of Mexican geoduck Panopea globosa (López-Landavery et al., 2021), whereas similar exposures delay metamorphosis in the Japanese geoduck P. japonica (Huo et al., 2019).
Transcriptome profiling of P. generosa has provided critical molecular insight on effects of low-pH exposure, highlighting molecular metabolic shifts over larvae-to-juvenile development and the informative capacity of transcriptomics for examining responses to low pH in this species (Timmins-Schiffman et al., 2020). There is also evidence that repeated exposures of juvenile P. generosa under OA conditions elicits compensatory growth and metabolism (Gurr et al., 2020) and differential DNA methylation (Putnam et al., 2017). In light of this, we posited that pre-exposure, or priming, generates changes in transcriptome profiles under repeated encounters. Specifically, we hypothesized that there will be frontloading of distinct gene functions and pathways underpinning adaptive phenotypes from early-life priming.
To this end, we examined gene expression patterns underlying high-pCO 2 priming in postlarval P. generosa that resulted in larger (tissue biomass and shell length) juvenile clams with reduced total antioxidant capacity after repeat exposures (Gurr et al., 2021). Gene expression data were analysed for juvenile P. generosa acclimated at the pediveliger stage under ambient and moderately elevated pCO 2 conditions for 110 days (day 0) and then were subsequently re-exposed in a reciprocal fashion to a second exposure of either ambient, moderately elevated pCO 2 , or severely elevated pCO 2 for 7 days (day 7), a 7-day return to ambient pCO 2 (day 14), and a third 7-day exposure to ambient or moderately elevated pCO 2 (day 21; Figure 1). Gurr et al. (2021) review details regarding the rationale (i.e., timing and magnitude of pCO 2 conditions) and limitations of this hormetic framework.

| pCO 2 exposure experiment and tissue sampling
Larval Pacific geoduck clams were reared from gametes at the Jamestown Point Whitney Shellfish Hatchery following standard industry practice (i.e., live-algal feed regime, larvae runts culled periodically, etc.). Once animals reached settlement competency (30 days post-fertilization), pediveliger larvae were exposed to ambient and elevated pCO 2 conditions (921 ± 41 and 2870 ± 65 μatm) for an initial 110-day conditioning period targeting the metamorphic transition from pediveliger to the burrowing juvenile stage (N = 4 trays treatment −1 and N = 1.5 × 10 4 pediveliger geoduck per tray). The timing for primary exposure has naturalistic relevance, as postlarval development represents a transition from a free-swimming stage to a sedentary life in the benthos (Goodwin and Pease, 1989), where bacterial carbon mineralization and low buffering capacity elevate pCO 2 and decrease calcium carbonate saturation (Cai et al., 2011).
Survivorship over the pediveliger-to-juvenile transition was ~30% regardless of pCO 2 condition (4-5 × 10 3 juvenile geoduck per tray; Gurr et al., 2021). Juveniles acclimated under ambient and elevated pCO 2 , hereafter referred to as naïve and pre-exposed clams, were divided at equal density into 36 replicate vessels (N = 120 animals per vessel, N = 6 vessels per treatment), and subjected to a secondary 7-day period under three pCO 2 conditions (ambient pCO 2 = 754 ± 15 μatm, moderate pCO 2 = 2750 ± 31 μatm, severe pCO 2 = 4940 ± 45 μatm) followed by 7 days of ambient recovery (896 ± 11 μatm) before replicates were split into 72 vessels (N = 6 vessels per treatment) for a 7day third exposure in two conditions (ambient pCO 2 = 967 ± 9 μatm and moderate pCO 2 = 3030 ± 23 μatm; Figure 1); the time to reach target treatments (i.e., ambient to moderately elevated pCO 2 ) occurred more rapidly (~3 h) than the return to ambient conditions from elevated pCO 2 levels (~ 6-8 hr). Note these pCO 2 values are higher than pCO 2 in the open ocean because they are designed to be relevant to the native range of Panopea generosa, as they correspond to measurements at local estuarine sites and sediment conditions where the clams live (e.g., >2400 μatm and Ω arag < 0.4 in Hood Canal, WA: Feely et al., 2010;JCP et al., 2014; Ω arag 0.4-0.6 in subsurface sediments: Green et al., 2009). Geoduck were fed a live mixed-algae diet ad libitum with a programmable dosing pump (Jebao DP-4), targeting 5 × 10 4 cells ml −1 in each vessel. Furthermore, marine bivalves can re-establish acid-base homeostasis 24-48 h after exposure to acidified seawater (Michaelidis et al., 2005;Spicer et al., 2007), and therefore we considered a span of 7 days as sufficient to infer a stable state during exposure to elevated pCO 2 and ambient seawater.
Additional details on geoduck rearing and experimental conditions are outlined in Gurr et al. (2021). As previously described in Gurr et al. (2021), the pre-exposed clams displayed a phenotype of reduced total antioxidant capacity and increased shell growth and tissue biomass under subsequent pCO 2 challenges, which provided evidence supporting the pediveliger-to-juvenile window for adaptive developmental plasticity. In this study, samples were sequenced at the same time points as physiological samples collected in Gurr et al. (2021)

| Gene expression
Individual whole juvenile geoduck samples were homogenized in DNA/RNA shield (1 ml) with 0.25 ml 0.5-mm glass beads by vortexing for 1 min. Total RNA was extracted from whole tissue homogenate using the Quick-DNA/RNA Kit (Zymo) according to the manufacturer's instructions. RNA was quantified using the Qubit RNA Broad Range Assay Kit with fluorometer (ThermoFisher) and quality was ascertained using a 4200 TapeStation System (Agilent Technologies) to visualize ribosomal bands, keeping in mind that the bands comigrate in geoduck clams into a single peak. RNA samples (10 ng μl −1 ) were used for TagSeq, a 3′ short transcript method that allows costeffective and accurate estimation of transcript abundances relative to traditional RNAseq (Lohman et al., 2016). Library preparation and sequencing of the 141 samples was conducted at the University of Texas Austin, Genomic Sequencing and Analysis Facility. Sequencing was completed on two lanes of an Illumina NovaSeq 6000 SR100, targeting standard coverage of 3-5 million 100-bp single-end reads per sample. Raw TagSeq reads were trimmed of Illumina adapters, poly-A sequence and quality filtered with fastp (Chen et al., 2018); quality control for filter optimization was completed using multiqc (Ewels et al., 2016). Following quality control, reads were mapped to the P. generosa reference genome (Putnam et al., 2022) using hisat2 F I G U R E 1 Experimental design for sampling. Line thickness indicates the pCO 2 treatments. Whole individual geoduck samples were snap frozen in liquid nitrogen at ~10:00 a.m. on each of days 0, 7, 14 and 21 and subsequently sequenced for gene expression (italicized numbers). White and grey geoduck cartoons correspond to primary exposure under ambient (naïve) and moderate elevated pCO 2 (pre-exposed), respectively. This figure is adapted from Gurr et al. (2021). (Kim et al., 2015) with a mapping efficiency of ~30%; the majority of unmapped reads aligned to ribosomal rRNAs (i.e., 18S and 28S) in geoduck clams (i.e., Panopea globosa). stringtie2 (Pertea et al., 2015) was used to quantify reads and assemble a count matrix (using prepDE. py) for analysis in R version 3.5.1 (https://www.r-proje ct.org). All code and data are publicly available (Gurr et al., 2022;doi:10.5281/ zenodo.6908630).

| Gene expression analysis
We chose a combined analytical approach to first test for gene expression frontloading (sensu Barshis et al., 2013) and fine-tuning or responsiveness to pCO 2 change by examining modules of coexpression via network analysis and their functional enrichment ( Figure 2) and second to test for differential expression in pairwise treatment groups at each time point (i.e., differentially expressed genes; DEGs) and their functional enrichment.
First, Weighted Gene Co-expression Network Analysis (wgcna; Zhang & Horvath, 2005;Langfelder & Horvath, 2008) was used to identify gene expression patterns and test for frontloading of transcripts (constitutive changes in expression, Figure 2c), or finetuning of transcripts under pCO 2 change (dynamic changes across repeated exposures), due to conditioning history by calculating gene co-expression modules at each time point. Co-expression network construction allows an assessment of broad expression-level directionality and the influence of compounding treatment history (e.g., initial × second exposure × third exposure), as opposed to typical approaches of pairwise differential expression analysis. Second, we used deseq2 to test for differential expression in pairwise group comparisons. While pairwise group comparisons are not able to examine the interactive effects, pairwise models investigated the effects of initial acclimation (ambient vs. moderate), second exposure (ambient vs. moderate, ambient vs. severe, and moderate vs. severe) and third exposure (ambient vs. moderate) and grouped contrasts to determine changes in gene expression due to cumulative (not interactive) pCO 2 treatment history.
Gene expression in response to pre-exposure and repeated encounters was analysed with wgcna using the bioconductor (version3.13) package "wgcna" (version 1.70-3; Langfelder & Horvath, 2008, Langfelder & Horvath, 2012 in R (R Core Team, 2021) to assess co-expression patterns (Zhang & Horvath, 2005). with the main effect of primary history only). To determine the functional patterns associated with pCO 2 treatment, co-expression modules significantly correlated with treatment variables were assessed for significantly enriched "molecular function" and "biological process" GO terms (p < .05) and GOslim assignment was used to place significantly enriched GO terms into hierarchical bins of broader function. Lastly, Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa & Goto, 2000) was applied to understand higher-level functional processes of co-expression modules and the web interface "keggmapper" (Kanehisa and Sato, 2020) was used to manually investigate genes involved in enriched pathways. Additional details regarding read filtration and gene expression analysis (wgcna and deseq2) are provided in the Supporting Information.
Co-expression module-treatment associations significantly correlated with the primary exposure period ( Figure 2a) were categorized as "frontloaded" or "fine-tuned" with respect to their expression patterns. Gene sets were identified as either consistently higher in expression through time (Figure 2b,c) were identified as putatively frontloaded, and those that were dynamic across repeated exposures ( Figure 2b) were identified as putatively fine-tuned. Quantitatively, gene frontloading was assessed as described in Barshis et al. (2013), where frontloaded genes are defined as those with a higher constitutive expression due to priming and less responsive to a subsequent environmental challenge. The term "frontloaded" was therefore assigned to those genes with higher expression under ambient pCO 2 by pre-exposed clams ("control ratio"; Figure 2c y-axis) and lower response ratio by pre-exposed clams ("response ratio" = expression in elevated pCO 2 : expression in ambient pCO 2 ; Figure 2c x-axis).
Frontloaded genes can then be defined as those with control ratio >1 and response ratio <1 (i.e., top left quadrant Figure 2c). A oneway ANOVA was used to test differences in vst-normalized gene expression of frontloaded genes between primary × second pCO 2 treatments. In contrast to frontloaded genes, fine-tuned expression was defined by first identifying the genes in modules correlated with exposure to the primary exposure period, and then second by tracking these genes underlying the GO terms that were uniquely respon-

| Co-expression network analysis overview
Network analysis using wgcna (Langfelder & Horvath, 2008;Zhang & Horvath, 2005) resulted in groups of genes that shared expression at a time point, termed modules. These modules are each named with a colour based on the convention of Horvath and Langfelder (Langfelder & Horvath, 2008;Zhang & Horvath, 2005). Given our goal to examine the interaction of repeated exposures, we analysed

| Day 0
Network analysis after the 110-day acclimation period (primary exposure, day 0) resulted in six co-expression modules with between 414 and 2818 genes ( Figure S1), excluding the module containing only one gene (the module "grey" which contains unassigned genes).
One module (Day 0 midnightblue) was significantly correlated with F I G U R E 2 Conceptual guide on our statistical approach to identifying gene frontloading and fine-tuning using co-expression analysis. White and grey geoduck cartoons represent the averaged expression data of naïve and pre-exposed clams. Shown is the directionality of module-treatment associations strongly correlated with the primary exposure period (a), containing groups of genes with similar expression patterns across time points. These putatively frontloaded genes were investigated for continuous or transient functional enrichment (b). Further, constitutive frontloading (sensu Barshis et al., 2013), defined by the quantitative criteria of genes with higher constitutive expression and lower response to reciprocal exposure, was used following co-expression network analysis. Genes that fall in the upper left quadrant are thus defined as "frontloaded" in their expression (c). [Colour figure can be viewed at wileyonlinelibrary.com] the primary exposure treatment and represented genes with higher vst-normalized gene expression values in the naïve clams ( Figure S1).
Note that the authors of wgcna recommend larger sample matrices (e.g., >15 samples) to yield robust and biologically relevant data (Langfelder & Horvath, 2008); because eight samples were sequenced on day 0, these results must be interpreted with caution and are presented in Figure S1.

| Day 7
Modules correlated with initial treatment Network analysis of expression following second exposure (day 7) resulted in five co-expression modules containing between 304 and 3951 genes ( Figure S2A). Of these five modules, two were significantly correlated with the primary exposure treatment (Day 7 module brown, and Day 7 module yellow Figure S2A). Day 7 module brown (N = 862 genes) showed higher vst-normalized gene expression values in the naïve clams relative to the pre-exposed clams, whereas Day 7 module yellow (N = 610 genes) had higher vstnormalized gene expression in the pre-exposed clams ( Figure S2B).

Modules correlated with interactions of treatments
One module was only significantly correlated with the combined primary and second treatments (Day 7 module green; Figure S2) and suggests that primary acclimation treatments and second exposures under moderate and severe pCO 2 lead to divergent expression patterns. Day 7 module green (N = 304 genes) showed low gene expression by pre-exposed clams re-exposed under moderate pCO 2 , whereas naïve clams had high expression values when exposed under severe pCO 2 ("MM" < "AS"; Figure S2). This module was enriched for cellular nitrogen compound metabolic process and included transcription factors and other regulators of transcription.
There were no pathways enriched for Day 7 module green.

Frontloaded genes
Of the full Day 7 module yellow (N = 610 genes), naïve clams increased expression of 346 and 184 genes under moderately elevated and severely elevated pCO 2 , of which 243 and 138 genes were assigned as "frontloaded" by pre-exposed clams; 106 genes were

| Day 14
Modules correlated with initial treatment Day 14 module black, representing genes highly expressed by pre-exposed clams, was enriched for "pentose phosphate pathway" (N = 5), "proteasome" (N = 6), "glycolysis/gluconeogenesis" (N = 7), "biosynthesis of amino acids" (N = 9) and "carbon metabolism" (N = 12; Table S1 and Figure 3). Gene families associated with F I G U R E 3 GO and pathway enrichment analysis of modules representing higher expression by naïve clams (white geoduck symbol, wgcna modules = Day 7 brown, Day 14 brown, and Day 21 blue and magenta) and pre-exposed clams (grey geoduck symbol, wgcna modules = Day 7 yellow, Day 14 black, and Day 21 yellow). Branches show enriched hierarchical GO terms ("bp" = biological process; "mf" = molecular function) unique to naïve and pre-exposed clams (solid and dashed lines) and those that were shared (wavy lines). GO terms are divided into two major groups: (i) "persistent terms" are enriched functions throughout the experiment; (ii) "transient terms" are enriched unique to or partially shared between the second exposure and ambient recovery period as "stress-induced and recovery," ambient recovery and third exposure period as "recovery and preparatory regulation,", and both elevated pCO 2 periods as "rapidly induced under stress exposures." Enriched pathways (p adj < .05) are shown for each module representing higher expression by naïve and pre-exposed clams.

Modules correlated with interactions of treatments
Four modules were significantly correlated with the combined primary and second treatments (Day 14 modules brown, black, pink and magenta; Figure S3A) and each suggests that primary acclimation treatments and second exposures under moderate and severe pCO 2 lead to divergent expression patterns. Of these four modules, two were not significantly correlated with the initial acclimation treatment (Day 14 module pink, and Day 14 module magenta; Figure S3A). Day 14 module magenta (N = 336 genes), representing genes with low expression values by pre-exposed clams re-exposed under severe pCO 2 , was enriched for "spliceosome" (N = 8; Table S2) involving alternative splicing proteins. GO analysis of Day 14 module magenta showed enriched terms in GOslim categories for cellular nitrogen compound metabolic process and ion binding ( Figure S8B) and included components of the spliceosome, E3 ubiquitin-protein ligases, tubulin proteins and Ras-related proteins. Day 14 module pink (N = 391 genes), representing genes highly expressed by preexposed clams re-exposed under severe pCO 2 ( Figure S3), did not show pathway enrichment (Table S1), but GO terms were enriched in GOslim categories of cellular nitrogen compound metabolic process, DNA binding, enzyme regulator activity and oxidoreductase activity ( Figure S8B), including homeodomain transcription factors, superoxide dismutase and Rho-GTPase activating proteins.

| Day 21
Modules correlated with initial treatment Network analysis of gene expression following the third exposure (Day 21) resulted in nine co-expression modules containing between 451 and 1753 genes ( Figure S4). Of these nine modules, three were significantly correlated with the primary exposure treatment (Day 21 modules blue, magenta and yellow; Figure S4A). Day 21 modules blue (N = 1537 genes) and magenta (N = 241 genes) showed higher vst-normalized gene expression values in naïve clams, whereas Day 21 module yellow (N = 926 genes) showed higher expression values in the pre-exposed clams ( Figure S4B).
Day 21 module blue, representing genes highly expressed by naïve clams, was enriched for "fatty acid degradation" (N = 11), "fatty acid metabolism" (N = 13), "endocytosis" (N = 31), "peroxisome" (N = 15) and "lysosome" (N = 25; Table S1 and Figure 3), involving the same genes and gene families enriched in earlier modules with the same expression pattern (i.e., Day 7 module brown and Day 14 module brown). Day 21 module magenta was enriched for "protein processing in endoplasmic reticulum" (N = 10; Table S2) and contained genes for protein quality control and trafficking and stress of the endoplasmic reticulum. All enriched GO terms for Day 21 modules blue and magenta were in GOslim categories that were also enriched in earlier modules representing the same expression pattern (i.e., Day 7 module brown and Day 14 module brown) with the exception of enzyme binding and immune system response ( Figures S5A, S6 and S9). These additional terms included genes for endosomal cargo trafficking and transport, regulation of ion channels, serine/threonine-protein kinases, E3 ubiquitin ligases and zinc finger proteins.
Day 21 module yellow, representing genes highly expressed by pre-exposed clams, was enriched for "mitophagy" (N = 7; Table S1 and Figure

Modules correlated with interactions of treatments
Six modules were significantly correlated with the combined primary, second and third treatments (Day 21 modules blue, magenta, red, black, pink, turquoise; Figure S4). Of these six modules, two were significantly correlated with the primary exposure treatment (Day 21 modules blue and yellow) and four were only correlated with the combined treatment history (Day 21 modules red, black, pink, turquoise; Figure S4A). Day 21 module black (N = 660), representing genes expressed at a low abundance by naïve animals under moderately elevated pCO 2 during the third exposure ( Figure S4A), was enriched for "ribosome biogenesis in eukaryotes" (N = 15), "spliceosome" (N = 19), "RNA transport" (N = 19), and "protein processing in endoplasmic reticulum" (N = 13; Table S1). There were no pathways enriched for Day 21 modules pink (N = 451 genes), red (N = 713 genes) or turquoise (N = 1753 genes; Table S1).

| Differential gene expression
Differential gene expression analysis on Day 0 resulted in 14 DEGs between treatments with higher expression level of 11 genes by naïve clams and three genes by pre-exposed clams (Table S2 and Subsequent exposures on days 7, 14 and 21 showed greater differential expression due to pCO 2 priming than second or third pCO 2 treatments (  Figure 5). In summary, the majority of DEGs resulting from primary exposure in this study were upregulated in the naïve phenotype (70%), especially following ambient recovery (85% of upregulated DEGs; Table S2 and DEGs (e.g., mucin-1, chitotriosidase-1 and defensin; Table S3). Persistent DEGs with higher expression level by naïve clams notably differed in their functional annotation for response to stress, lipid and calcium binding, and protease inhibition (apolipoprotein D, regucalcin, mucin-1, four-domain protease inhibition, etc.; Table S3). Persistent DEGs with higher expression level by pre-exposed clams, although fewer, were additionally associated with antimicrobial activity, cobalt transport (cobalamin) and protease inhibition (defensin, CD109 antigen, and BPTI/Kunitz domain-containing protein 4-like;

| DISCUSS ION
Post-larval acclimation to elevated pCO 2 affected transcriptome profiles. Specifically, priming to moderately elevated pCO 2 involved F I G U R E 5 Effect of ambient and moderately elevated pCO 2 priming on differential gene expression. Venn diagrams show the number of genes upregulated (a, higher expression level in naïve clams) and downregulated (b, higher expression level in pre-exposed clams) and unique to and shared between experiment days; featured counts in bold represent persistent DEGs shared among days 7, 14 and 21. Segment charts represent significant "molecular function" GO terms (−log 10 p-value) in greyscale for time points.

| Naïve profile: endogenous lipids supply high transcriptional demand
A growing body of research suggests that environmental changes, such as elevated pCO 2 /low pH conditions, increase energy partitioning toward protein biosynthesis (Langenbuch & Pörtner, 2002;Pan et al., 2015), conferring energetic tradeoffs for somatic growth and storage retention (Sokolova, 2013;Stumpp et al., 2011). Our results support this, as a key difference between transcriptome profiles was a higher transcript abundance attributed to fatty acid degradation and glutathione components in the naïve phenotype (Figures 3   and 5). Persistently enriched pathways included peroxisome activity (β-oxidation), acetyltransferase to mitochondria (carnitines) and bioremediation of free radicals illustrating elevated use of endogenous metabolic fuel to satisfy greater transcriptional demand.
Mobilization of endogenous reserves, primarily lipids, is essential to meet energetic requirements of early development Waldbusser et al., 2013), but also plays a vital role in rapid energy provisioning during stress exposure  and may be advantageous to sustain acid-base homeostasis (Evans et al., 2017). Naïve juveniles were enriched for carbohydrate metabolism (i.e., pentose phosphate pathway, glycolysis and carbon metabolism) prior to pCO 2 challenge, further representing greater energetic requirements and the reliance for lipid sources to satisfy demand under pCO 2 change.
Across multiple marine calcifiers, exposure to elevated pCO 2 causes lipid loss and altered fatty acid metabolism coupled with shell malformations and delayed settlement competency (Dickinson et al., 2012;Liu et al., 2020;Talmage & Gobler, 2010;Timmins-Schiffman et al., 2020). For instance, elevated pCO 2 affects shell biomineralization and fatty acid metabolism in the pearl oyster Pinctada fucata  and reorganizes the lipid profile in purple-hinge rock scallop Crassadoma gigantea (Alma et al., 2020). Gurr et al. (2021) found naïve geoduck clams decreased organic biomass and shell length under elevated pCO 2 , linking physiological responses with gene-expression patterns of naïve clams in this study. Similarly, marine invertebrates are known to increase lipid degradation and peroxisome activity under elevated pCO 2 , such as higher expression levels of long-chain specific acyl-CoA dehydrogenase in the coral Acropora millepora and barnacle Balanus amphitrite (Wong et al., 2011;Kaniewska et al., 2012) and peroxiredoxins and carnitine O-acetyltransferase in larval oysters Crassostrea virginica and Crassostrea hongkongensis (Tomanek et al., 2011;Dineshram et al., 2015). In contrast, elevated pCO 2 may not affect fatty acid metabolism (Matson et al., 2012;Timmins-Schiffman et al., 2014), or may interact with multiple stressors on lipid use (e.g., dietary restriction; Gibbs et al., 2021), which is a testament to an array of contingencies affecting metabolic shifts (e.g., species, timing, stress type[s] and intensity). Analysis of the lipidome (totality of lipids in an organism) can assess the importance of lipid catabolism on physiological success (Laudicella et al., 2020) and future efforts should consider the tissue-specificity of proteomic and gene expression patterns (e.g., Elowitz, 2002;Wei et al., 2015), requiring fine-scale sampling in contrast to whole-tissue homogenates sequenced herein.
Altogether, the gene expression patterns of naïve Panopea generosa suggest that fatty acid degradation is a primary response to OA to ensure short-term survival and may compensate for the additional transcriptional requirements; however, depletion of endogenous storages confers an unsustainable mismatch between energy demand and supply if the pCO 2 challenge exacerbated or persisted (i.e., "pessimum" range; Sokolova et al., 2012;Sokolova, 2021). Thus, the transcriptomic profile of pre-exposed P. generosa illustrates the importance of early-life environmental cues for eliciting molecular resistance. Beyond the scope of this study, standing and cryptic genetic variation may underlie adaptive capacity and heritable plasticity to environmental change (Bitter et al., 2019). For example, survival of the Yesso scallop Patinopecten yessoensis (Yang et al., 2021) and normal development of the purple sea urchin Strongylocentrotus purpuratus (Pespeni et al., 2013) under elevated pCO 2 is linked to allele variation of candidate genes associated with energy and lipid metabolism, respectively. Future studies should further examine genomic markers affecting transcriptome profiles and selection.
A substantial increase in DEGs during repeat exposures, relative to immediate post-acclimation, suggests that phenotypic differences are most evident upon environmental change. In particular, a greater magnitude of change in gene expression occurred during ambient recovery ( Figure 5 and Table S2). Similarly, mussels

| Primed profile: fine-tuned response to intermittent OA
Pre-exposed P. generosa expressed fewer genes at a higher abundance relative to naïve geoducks, albeit functional annotation of these genes showed diverse regulatory processes in cellular quality control and homeostasis (Figure 3). A general decrease in transcription may be attributed to adaptive benefits under environmental change (Bultelle et al., 2021). For instance, mussels Mytilus californianus decrease gene expression when acclimated to dynamic thermal stress, as opposed to acute isothermal conditions, highlighting a lower transcriptional demand during intermittent exposures (Connor & Gracey, 2020).
generosa also expressed higher levels of genes essential for protein turnover, such as 26S proteasome, E3 ubiquitin ligases and caspase (Goldberg, 2003;Voges et al., 1999), that are otherwise unaffected by low pH in other species (i.e., C. virginica and Mercenaria mercenaria; F I G U R E 6 Summary of transcriptome patterns expressed by pre-exposed and naïve juvenile geoduck clams under repeat exposure to elevated pCO 2 . Putative mechanisms are shown based on synthesis of the enriched terms (functions and pathways) and differential expression in this study, alongside phenotypic data described in Gurr et al. (2021) (indicated by asterisks and italics; "TAOC" total antioxidant capacity); recommended topics for future research are proposed ("?"). Götze et al., 2014) probably due to energetic limitations of exposure to environmental change (Ivanina et al., 2016). Pre-exposed P. generosa increased expression of genes for energy metabolism (NADH dehydrogenase, cytochrome c oxidase and ATPase) and biosynthesis (pentose phosphate pathway) during ambient recovery, suggesting an opportunistic or pre-emptive increase in energy production under optimal conditions (Figures 3 and 6).
Stimulation of the electron transport chain increases energy production under low-pH conditions (Evans et al., 2017), but altered expression of mitochondrial complexes may also confer metabolic suppression. For example, under elevated pCO 2 , geoduck P.
globosa (López-Landavery et al., 2021), Pacific oyster C. gigas (Wei et al., 2015) and eastern oyster C. virginia (Chapman et al., 2011) upregulate NADH dehydrogenase, suggesting increased ATP production, whereas lower levels of gene expression for cytochrome c oxidase in C. gigas (Dineshram et al., 2012) and ATP synthase in C. hongkongensis (Dineshram et al., 2013) suggest metabolic suppression. Thus, increased expression of mitochondrial complexes by pre-exposed geoduck clams illustrates an opportunistic increase in energy production. Moreover, enrichment for glycolysis and the nonoxidative pentose phosphate pathway suggests the pre-exposed clams metabolized carbohydrate fuels and increased nucleotide biosynthesis during ambient recovery. This demonstrates expeditious rescue by primed clams, as naïve clams were enriched for carbohydrate metabolism under ambient priming but not post-challenge during ambient common garden conditions. GO term enrichment during ambient recovery also included iron binding proteins and phenoloxidases ("oxidoreductase activity"; Figure S6). Expression of ferric-chelate reductase may improve iron homeostasis and prevent excess iron-induced toxicity (Li et al., 2019), converting ferric iron to an "active" electron-donor state (ferrous iron) required for biological processes (Connolly et al., 2003). Since antioxidants were not abundantly expressed by pre-exposed clams, excess iron-induced toxicity (Fenton reaction enhancing free radicals) was probably negligible.
In contrast, stressed mud snails Littorea littoriea increase expression levels of antioxidants and ferritin (Larade & Storey, 2004) essential for storing iron, suggesting taxon-specific flux of iron constituents during stress. Lastly, laccase and tyrosinase are phenoloxidases of growing interest as biomarkers of immune response and detoxification (Luna-Acosta et al., 2017) and were expressed during the ambient recovery period by pre-exposed clams. Future study is needed to determine the role of divergent transcriptome profiles during intermittent acidosis in marine invertebrates.
Lastly, immunomodulation and signalling were key transcriptional components of pre-exposed P. generosa following the acclimatory period, including genes for activation and inhibition of NF-kappa β (Figures S5B and S6; e.g., mitogen-activated protein kinase, toll-like receptors 2, 3 and 4, TNF receptor-associated factor 6, MyD88, B-cell lymphoma 3 protein, and death-associated inhibitor of apoptosis 2), a transcription factor involved in immune deficiency signalling cascade in defence of pathogens (Leulier et al., 2006).
Similarly, immune priming under abiotic (i.e., temperature) and biotic (i.e., synthetic and nonsynthetic viral exposure) challenges bolsters antiviral defences associated with upregulated expression of apoptotic, signalling and immune-related genes in the Pacific oyster C.
gigas (Delisle et al., 2020;Lafont et al., 2020) and Pacific abalone Haliotis discus hannai (Zhang et al., 2022). A growing body of research highlights the general importance of NF-kappa β for the innate immune response in bivalves Li et al., 2015) and elevated pCO 2 can have synergistic and antagonistic effects on immunomodulation (Cao et al., 2018;Castillo et al., 2017). For example, increased levels of NF-kappa β in the mussel Mytilus coruscus may improve immune defences to compensate for weakened shell strength under low pH (Zhao et al., 2020). Moreover, the blood clam Tegillarca granosa decreases gene expression of NF-kappa β signalling during hypercapnia, rendering clams more susceptible to disease . After an initial stress encounter, Mytilus galloprovincialis reduces transcription of immune-related proteins, but insufficient to counteract decreased growth (Détrée & Gallardo-Escárate, 2018).
Altogether, early-life experience heightened critical signalling and immune-related proteins potentially enhancing resilience to subsequent pCO 2 change in primed clams. Considering that pre-exposed P. generosa grew larger than naïve clams (Gurr et al., 2021), earlylife priming to elevated pCO 2 probably affects organismal responses to subsequent encounters demonstrated herein from "fine-tuned" transcript abundance during repeated exposure.

| Transcriptional control suggests "memory" post-acclimation
Understanding how the environment triggers biological responses that lead to gene expression regulation is key to understand- Stress "memory," or stored/imprinted/regulatory information from initial stress enhancing robustness to future encounters, has largely been studied as a plant-based phenomenon (Bruce et al., 2007), with growing support in invertebrate models. Molecular mechanisms underpinning memory may manifest as nongenetic markers, transcription factors and key signalling metabolites with cascading implications for performance. For example, sustained expression of the transcription factor Nrf2 co-occurs with improved antioxidant defence systems in cold-primed tunicates Ciona robusta (Li et al., 2020). In our study, pre-exposed P. generosa frontloaded several transcription factors, suggesting an influence of stress history on general regulatory transcription. Moreover, histone methyltransferases (HMTs) were expressed at a higher abundance by pre-exposed clams (i.e., SETD5, ASH1L and NSD2) each affecting histone H3 tri/dimethylateion of lysine residue 36 (H3K36me3 and H3K36me2; Greer & Shi, 2012), a chromatin-carrying marker for recruitment of gene-body DNA methylation (Dhayalan et al., 2010;Nanty et al., 2011), alternative splicing (de Almeida et al., 2011), and coparticipating in histone acetylation (Osipovich et al., 2016). Gene expression frontloading can enhance the ability for marine invertebrates to cope with dynamic environments (Barshis et al., 2013;Dong et al., 2008;Shiel et al., 2017), and may be rooted in epigenetic processes. For example, normal development of larval oysters C. honkongensis under low pH co-occurred with upregulation of HMTs and fine-tuned transcription (Dineshram et al., 2015). Thus, the continuous expression of histone modifiers by pre-exposed P. generosa in this study suggests a rapid and inducible mechanism to support gene frontloading. Altogether, sustained/accumulated HMTs may regulate DNA accessibility for transcription and contribute to the fine-tuned and responsive gene expression and emergent phenotype in juvenile geoduck clams, consistent with the well-established role of histone modifications in stress "memory" and improved performance (Mozgova et al., 2019).
Genome-wide epigenetic and post-translational modifications can mediate phenotypic variation (Anastasiadi et al., 2017;Liew et al., 2020;Putnam et al., 2016) and fine-tune transcription (Liew et al., 2018), although in some cases this mechanism is subtle (Downey-Wall et al., 2020). P. generosa with a history of low pH exposure have demonstrated epigenetic signatures of differentially methylated genes linked to a beneficial phenotype of compensatory growth and resilience when challenged with low pH again, in comparison to more sensitive clams exposed to low pH for the first time (Putnam et al., 2017). As an expansion of this finding, expression of HMTs may control accessibility of DNA with cascading effects on essential biological processes (e.g., signal transduction, cell proliferation, growth and cell death; Greer & Shi, 2012). Moreover, pediveliger larval oysters C. hongkongensis exposed to OA show hypermethylation of genes related to lipid metabolism in contrast to hypomethylation of genes related to carbohydrate metabolism (e.g., NADH dehydrogenase and ATP synthases), suggesting that lipid fuel acts as an alternative energy source to cope with elevated pCO 2 (Lim et al., 2021). Further research is needed to elucidate molecular mechanisms underpinning gene-expression regulation, and expanded efforts require long-term tracking and interdisciplinary approaches (multi-omics) to understand how plasticity affects species and populations with different susceptibilities (Fox et al., 2019).
Moreover, standing genetic variation may drive adaptive directional selection under environmental change, as evidenced by sufficient genetic variation for rapid evolution in the marine molluscs (Bitter et al., 2019;Thomsen et al., 2017;Yang et al., 2021) and echinoderms (Pespeni et al., 2013) to cope with elevated pCO 2 .

| CON CLUS ION
In this study, we investigated the transcriptome profiles of juvenile geoduck clams post-acclimation and under intermittent pCO 2 change to test the hypothesis that priming during a developmental window elicits transcriptional frontloading and an expeditious response to environmental challenges. In the absence of priming, naïve clams exhibited higher transcriptional demand during experimental exposures attributed to a putative catabolic shift favouring primarily lipids. In contrast, acclimatization conferred gene frontloading and gene-expression control, such that pre-exposed P. generosa contained lower levels of expression, but frontloaded and fine-tuned under subsequent pCO 2 /pH changes. The pre-exposed transcriptome profile was attributed to heightened expression for cellular quality control, signalling, transcriptional modifiers, and stress and innate-immune response genes under subsequent pCO 2 challenge and putative markers of enhanced energy production during recovery. Altogether, this study demonstrates the importance of geneexpression regulation on positive developmental acclimatization.

ACK N OWLED G EM ENTS
We thank the Jamestown S'Klallam Tribe and Kurt Grinnell for providing the animals and facilities for this research. We also thank management staff and technicians at the Jamestown Point Whitney

CO N FLI C T O F I NTE R E S T
The authors have declared no conflict of interest for this article.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5281/ zenodo.6908630.