Adult purple sea urchins (Strongylocentrotus purpuratus, N = 83) were collected on September 20, 2017 off of Goleta Beach, CA, United States (34° 24.840′ N, 119° 49.742′ W), at a subtidal site near the University of California Santa Barbara (UCSB) under California Scientific Collection permit SC-1223. Urchins were transported to UCSB where they were held in flow-through seawater tables at ambient seawater conditions for 6 days. After this initial acclimation period, urchins were randomly assigned to one of two acclimation treatments, upwelling (U: 1390 ± 307 μatm pCO₂ and 12.7 ± 0.5°C) or non-upwelling (N: 631 ± 106 μatm pCO₂ and 16.8 ± 0.2°C) (Table 1 and Supplementary Figure 1). These conditions represent end points of those experienced in their natural environment (Kapsenberg and Hofmann, 2016; Rivest et al., 2016). Each treatment was replicated across three 20-gal flow-through aquarium tanks with each holding 12 urchins. Conditioning of adults was maintained for just over 4 months and animals were fed locally collected kelp (Macrocystis pyrifera) once a week.

Adult acclimation conditions were controlled by a flow-through CO₂-mixing system (modified from Fangue et al., 2010) and Delta Star heat pumps regulated by Nema 4x digital temperature controllers (AquaLogic San Diego, CA, United States). Each treatment tank was fed water from two 20-gal glass reservoir tanks supplied with 0.35 μm-filtered, UV sterilized seawater (FSW) with a flow rate of 10 L/h (20 L/h total into each tank). Target pCO₂ levels were maintained in the reservoir tanks using mass-flow controllers (Sierra Instruments, United States) that mixed pure CO₂ with dry air that was filtered and CO₂ scrubbed.

Seawater chemistry measurements in each tank were taken every 1–2 days throughout the adult acclimation period (Table 1 and Supplementary Figure 2). Temperature and salinity were measured daily using an Omega HH81A thermocouple and YSI-3100 conductivity meter, respectively. pH was measured using a UV spectrophotometer (Shimadzu UV-1800) and m-cresol purple indicator dye (Sigma Aldrich) according to standard operating procedure (SOP) 6b (Dickson et al., 2007). Total Alkalinity (TA) was measured by titration following SOP 3b (Dickson et al., 2007). Water samples for total alkalinity were preserved prior to titration with 0.02% saturated mercuric chloride. CO₂ calc (Robbins et al., 2010) was used to calculate the carbonate chemistry parameters pCO2, Ωara, and Ωcal with the equilibrium constants from Hansson (1973) refit by Dickson and Millero (1987; Table 1).

On January 30th, 2018 spawning of the acclimated adult urchins was induced by injection of 0.53 M KCl into the coelom. Sperm was collected by pipetting from the area around the gonopores of male urchins kept out of water and stored dry in 1.5 mL tubes on ice until activation. Eggs were collected from each female by inverting the female over a beaker filled with FSW so that the gonopores were submerged. Eggs were checked visually for quality under a microscope and counted. Females with eggs of abnormal shape or color or significant numbers of eggs containing visible germinal vesicles were excluded. Sperm from potential males was activated in FSW, checked for motility, and then used to fertilize a subset of eggs from each female to ensure male-female compatibility. Eggs from 9 females per treatment (N or U) were pooled in equal numbers and fertilized with the sperm from one high quality non-upwelling conditioned male (Supplementary Figure 1). As this study focused solely on maternal contributions to intergenerational plasticity, one male was used to sire all offspring to minimize the genetic diversity that might generate noise in the analysis. As such, all larvae in the experiment were a mix of full and half siblings. Dilute, activated sperm was slowly added to each pool of eggs until at least 95% fertilization was reached. The embryos from each treatment pool were divided into 3 non-upwelling culture vessels (N: 537 ± 28 μatm and 17.2 ± 0.3°C) and 3 upwelling culture vessels (U: 1145 ± 87 μatm and 13.3 ± 0.3°C) (Table 1 and Supplementary Figures 1, 2) in a flow-through seawater system for a final concentration of 10 embryos/mL (120,000 embryos/vessel). Therefore, embryos were raised either in the same treatment as the maternal acclimation or the reciprocal condition. This yielded four treatment groups, each with 3 replicate culture vessels: embryos from mothers conditioned to non-upwelling raised under non-upwelling (NN), embryos from mothers conditioned to non-upwelling raised under upwelling (NU), embryos from mothers conditioned to upwelling raised under non-upwelling (UN), embryos from mothers conditioned to upwelling raised under upwelling (UU).

Temperature and pCO₂ conditions under which the embryos developed were maintained with the same CO₂ mixing system and heat pumps used in the adult acclimation. Two reservoir tanks fed 6 culture vessels per treatment (non-upwelling or upwelling developmental conditions) at a rate of 6 L/hr using irrigation drippers (DIG Corporation). Each culture vessel consisted of two nested buckets with a 12L capacity in which the inner bucket had holes covered with 30 μm mesh, allowing water to flow from the inner to the outer bucket and overflow without losing any embryos. Flow-through and mixing was aided in each bucket by an acrylic paddle attached to a 12V motor. Following the same methods outlined for the adult acclimation treatments, temperature and pH were measured daily for each culture vessel to confirm uniform conditions across each treatment. Salinity and TA were measured daily as well (Table 1 and Supplementary Figure 2). Embryos were cultured in their respective treatments until the prism stage, an early echinopluteus larval stage, at which point they were sampled to assess morphometrics, DNA methylation, and transcription.

The prism stage was defined by the archenteron merging to one side of the body and becoming tripartite, the first development of skeletal rods, and a pyramid-like body shape. Prism larvae were preserved for morphometric analysis by adding 4% formalin buffered with 100mM NaBO₃ (pH 8.7) to an equal volume of larvae and FSW for a final concentration of 2% formalin. Samples of 1000 larvae were stored at 4°C until measured. Larvae were photographed on a compound microscope (Olympus BX50) with an attached digital camera (Motic 10 MP) using the Motic Images Plus software. Prism larvae were photographed from a lateral view where both the tip of the body rod and postoral rod were in focus. Photographs were calibrated to a stage micrometer. Spicule length was measured as the length from the tip of the body rod to the branching point of the postoral rod (N = 35 larvae/culture vessel) (Supplementary Figure 4). The effect of each fixed factor (maternal condition or developmental condition) on spicule length was examined using a 2-way ANOVA.

Genomic DNA was extracted using an established CTAB protocol involving phenol chloroform purification. One sample per culture vessel was prepared for sequencing to generate 12 total samples (3 per treatment × 4 treatments). Each sample contained a pool of 6,000 larvae. Therefore, for each treatment, three replicates of 6,000 larvae were extracted and subsequently sequenced, totaling 18,000 larvae per treatment. In the case of the UU treatment, a missing sample from UU2 led us to process two unique pools (6,000 larvae each) of the UU3 culture. Integrity and quality of genomic DNA was assessed using gel electrophoresis and quantified using a Qubit fluorometer 3.0 (Life Technologies). Genomic DNA (∼2 μg per sample) was sent to the UC Davis Genome Center for RRBS library prep using the Diagenode Premium RRBS kit and sequencing on the Illumina 4000 with 50 bp single end reads across one lane of sequencing. Raw fastq sequences were trimmed and filtered using TrimGalore (V0.5.0)¹, which utilizes cutadapt (V1.9.1) (Martin, 2011), specifying the -rrrbs option that trims one additional bp from the 3′ end of the read as is recommended for RRBS data. This trimming package removes low quality base pairs from the 3′ end of the read (Phred score = 20), Illumina adaptors and short reads. Read quality was visualized using FastQC. Bisulfite treated trimmed reads were mapped to a bisulfite converted genome [Strongylocentrotus purpuratus genome (V3.1)] using Bismark (V0.19.1) (Krueger and Andrews, 2011). Methylation output files were produced by running the bismark_methylation_extractor command and summary files were produced by the bismark2summary command.

Bismark coverage files were used for the subsequent analysis in the R environment (V3.4.2). Statistical analysis of CpG data was run using the package methylKit (Akalin et al., 2012). Samples were filtered by read coverage where base pairs with less than 10X coverage were removed as well as base pairs in the highest percentile (99.9). CpG sites were merged across samples into a methylBase object where base-pair locations were the same across all samples using the unite function. Two independent tests were run to identify differential methylation at each CpG site (DMCpG), one to determine DMCpGs by maternal condition (upwelling vs. non-upwelling) and a second by developmental condition (upwelling vs. non-upwelling). DMCpGs were identified using a logistic regression and significance was determined using p-values adjusted to q-values using a sliding linear model method (SLIM) (Akalin et al., 2012). A CpG was considered significantly differentially methylated (DMCpG) if the q-value < 0.01 and the difference in percent methylation was >25%. Multivariate analyses were performed using percent methylation values for each CpG site. Principal coordinates analysis (PCoA) using Manhattan distances was performed using the packages adegenet (Jombart, 2008) and vegan (Oksanen et al., 2013). To determine the proportion of variance explained by fixed factors a permutational multivariate ANOVA was run using the function adonis.

CpG annotation to genome features was performed using a custom R script. A CpG site was considered methylated if its median value across all samples was >1 (the site was methylated in over half of the samples). To examine percent methylation values per gene, genes with sequenced CpGs were filtered for those with at least 5 representative CpG sites within the gene region. The remaining genes had an average of 10.3 total CpGs represented within the gene region and percent methylation was calculated for each gene (methylated CpG sites/total CpG sites). Gene Ontology (GO) functional enrichment was performed using a Fisher’s Exact test on genes that contained at least one DMCpG for either experimental treatment (maternal or developmental) compared to all genes using the GO_MWU package (Wright et al., 2015).

Total RNA was extracted using TRIzol reagent (Invitrogen). One sample per culture vessel was prepared for sequencing to generate 12 total samples (3 per treatment X 4 treatments). Each sample contained a pool of 6,000 larvae. RNA quality and quantity were assessed using the Agilent TapeStation and Qubit fluorometer 3.0 (Life Technologies), respectively, and all samples had RIN values >8.0. Total RNA (∼2 μg per sample) was submitted to the UC Davis Genome Center for poly A enriched, strand specific RNA library preparation and sequencing on the Illumina 4000 with 50bp single end reads across one lane of sequencing. Raw sequences were filtered for quality and adaptor sequences were discarded using TimGalore (V0.5.0)¹ which utilizes cutadapt (V1.9.1) (Martin, 2011). Trimmed reads were mapped to the Strongylocentrotus purpuratus genome (V3.1) using hisat2 (V2.0.0) (Kim et al., 2015). Gene counts were determined using featureCounts, a component of the subread package (V1.6.2) (Liao et al., 2014).

All statistical analyses were performed in R (V3.4.2). The package arrayQualityMetrics was used to identify potential sample outliers (Kauffmann et al., 2009). Gene normalization and models for differential gene expression analysis were performed using DESeq2 (Love et al., 2014). A Wald test [design = formula(∼ treat_maternal + treat_dev)] followed by pair-wise contrasts was run to identify differentially expressed genes (DEGs) associated with either maternal condition (upwelling vs. non-upwelling) or developmental condition (upwelling vs. non-upwelling). Multivariate analyses were performed using regularized log transformed gene count data. Principal coordinates analysis (PCoA) using Manhattan distances was performed using the packages adegenet (Jombart, 2008) and vegan (Oksanen et al., 2013). To determine the proportion of variance explained by fixed factors a permutational multivariate ANOVA was run using the function adonis.

Using only highly expressed genes (15,242 genes, mean counts >10), regularized log transformed counts data were used for a weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008). WGCNA analysis clusters genes in a manner that is un-biased toward the experimental design. A similarity matrix of gene expression was calculated using Pearson correlations for all gene pairs across samples and a “signed” network was specified, which allows the direction of the expression change to be retained. Using a power adjacency function and a soft thresholding power of 9, connectivities were transformed from the expression correlations. Genes were hierarchically clustered based on topological overlap which identifies groups of genes, or modules, whose expression covaried across samples. The minimum module size was set to 30 and similar modules were merged when their module eigengenes were highly correlated (R > 0.80). Experimental conditions and physiological larval trait data were then correlated to expression modules (Supplementary Figure 5). Gene Ontology (GO) enrichment analysis was performed for each module using a Fisher’s Exact Test, by presence or absence of each gene within a module using the GO_MWU package (Wright et al., 2015).

In order to examine how plasticity in the methylome and transcriptome correlate to a plastic, organismal trait, we measured the effect of different conditions on spicule length, a trait associated with larval fitness that is known to be sensitive to high pCO₂ but not to temperature differences within the range manipulated in this experiment (Byrne et al., 2013; Padilla-Gamiño et al., 2013; Wong et al., 2018). For early life-history stages of echinoderms, the larval skeleton begins to form just following gastrulation (Okazaki, 1975) and is contingent on sufficient Ωara, and Ωcal to support biogenic calcification (Byrne et al., 2013). We found strong evidence that spicule length varied between development treatments (p_ANOVA < 0.001), but no evidence for intergenerational plasticity (Supplementary Figure 4). Specifically, developing under upwelling conditions (higher pCO₂ and lower temperature) reduced average spicule length by 30.8% relative to those that developed under non-upwelling conditions. Thus, in the case of this organismal trait, we observed that abiotic conditions experienced during development were more influential than the environmental history of the mothers during gametogenesis.

To identify the potential for ecologically relevant abiotic conditions to induce changes in DNA methylation, either in an inter- or intragenerational plasticity context, we characterized the larval methylome using RRBS sequencing. DNA methylation sequencing yielded an average of 24,800,247 reads per sample across 12 samples (Supplementary Table 1), with an average bisulfite conversion rate of 98.92%. After adaptor trimming, read mapping to the bisulfite converted genome resulted in an average of 38.13% uniquely mapped reads and 23.44% ambiguously aligned reads (Supplementary Table 1). Bismark revealed 164,062,595 ± 9,335,444 SE total cytosine calls across all sequenced samples. Of these total calls, three different types of cytosine methylation were quantified, CpG (32,803,060 ± 1,885,036 SE), CpH (26,641,240 ± 1,496,753 SE), and CHH (104,619,294 ± 5,954,888 SE). Of these methylation calls, an average of 22.12% of CpGs were methylated, with 1.72% and 1.41% of CpHs and CHHs methylated, respectively (Supplementary Table 1). After filtering low coverage and high coverage CpG sites, read coverage per CpG site was 30.26 ± 0.962 SE and percent methylation per base was 22.95% ± 0.3957 SE (Supplementary Table 1). After filtering and sample merging, 245,343 CpG sites were shared across all samples, which were used for subsequent analysis.

Using this set of shared CpG sites, we sought to identify whether there were significant changes in the methylome associated with maternal or developmental conditioning. Logistic regression models revealed multiple differentially methylated CpGs. In the context of maternal effects, models found 684 DMCpGs total, with 328 hypermethylated and 356 hypomethylated sites in the methylome when mothers were conditioned in upwelling conditions compared to non-upwelling (Supplementary Figure 3). With regard to developmental conditioning, models found 216 DMCpGs total, with 175 hypermethylated and 41 hypomethylated sites when larvae developed in upwelling conditions compared to non-upwelling conditions (Supplementary Figure 3). Only 13 of these DMCpGs were significant for both maternal and developmental conditioning of which the majority were hypermethylated with regards to developmental conditioning and hypomethylated with regards to maternal conditioning (Supplementary Figure 3).

Principal coordinates analysis revealed that variation in maternal conditioning separated the samples along PCo1 while developmental conditioning separated samples along PCo2 (Figure 1A). A permutational multivariate ANOVA revealed that 12.8% (p < 0.001) of the variance in percent CpG methylation was explained by maternal condition while 9.7% (not significant) of the variance was explained by developmental condition (Figure 1C). The interaction between these two factors explained 8.5% of the variance, but was non-significant (Figure 1C). However, the majority of variance was not explained by any of the factors in the model, an outcome that is similar to methylation work in corals (Dixon et al., 2018; Baums et al., 2019), suggesting this remaining epigenetic variance may be driven by genetics or spontaneous epimutations. In summary, we found that maternal conditioning, the abiotic conditions that female urchins experienced during gametogenesis, resulted in more pronounced changes in the methylome than when embryonic development occurred in different abiotic conditions.

In order to examine relationships between DMCpGs, transcription and environmental conditioning, we characterized gene bodies containing DMCpG sites. Research on DNA methylation in invertebrates has highlighted patterns between gene body methylation (GBM) and broad transcriptional patterns (Gavery and Roberts, 2013; Dixon et al., 2014, 2018). Therefore, we sought to identify genes with DMCpGs in comparisons of both maternal and developmental conditions. In addition, we performed functional enrichment tests, Gene Ontology (GO), to see if broad functional categories were over-represented in those genes. A total of 95,197 CpG sites were located within gene regions across 9,219 total genes and are further designated as gene specific CpGs (GS-CpGs). Across GS-CpG sites there was a strongly bimodal pattern in the percentage to which each site was methylated (Supplementary Figure 6A). For each GS-CpG site, we considered the site methylated if the site had methylated reads in over half of the samples (median > 1), a threshold that resulted in 41,940 methylated GS-CpG sites. Genes were filtered for only those that had at least five GS-CpG sites sequenced per gene, resulting in a total of 5,024 genes. Percentage methylation per gene also showed a bimodal pattern, with the majority of represented genes being 100% methylated (Supplementary Figure 6B).

We found a small fraction of total genes that contained DMCpGs. A total of 136 genes contained at least 1 DMCpG for maternal condition while 49 genes contained at least 1 DMCpG for developmental condition. Overall, a total of 258 DMCpGs were also GS-CpGs across 177 unique genes (Supplementary Table 2). Notably, the majority of these 177 genes include those that function in protein modification including hydrolase activity (Sp-Vps35, Sp-Abhd3, Sp-Abhd3, Sp-Asah2L3), protein phosphorylation (Sp-Dusp6/7/9, Sp-Mekk4a, Sp-Melkl, Sp-Frap_1, Sp-Ck1g), and ubiquitination (Sp-Ube3cL, Sp-Uba1, Sp-Traf3). In addition to genes with protein modification processes, there were also genes involved in binding (Sp-Myef2, Sp-Z371, Sp-Z163, Sp-Ciz, Sp-Traf3, Sp-Pcbp2) and modification of nucleic acids (Sp-Endrvt33, Sp-Tsen34L_1). Cytoskeletal components such as myosin and kinesin are also represented across all genes containing DMCpGs (Sp-Non/MmyIIhcph, Sp-Unk_109, Sp-Krp85_1, Sp-Kif13A2, Sp-Kif1a, Sp-Kif13BL) as well as DNA binding proteins such as transcription factors (Sp-FoxK, Sp-Mlxip, and Sp-Myc) (Supplementary Table 2). GO enrichment tests were run on genes with DMCpGs for maternal condition, developmental condition, and all genes containing a DMCpG. Only tests for maternal DMCpG genes and all genes yielded significant GO enrichments, with the majority of the terms being the same (Supplementary Table 3). Among all genes with DMCpGs, there was significant GO enrichment of terms that reflect functions mentioned above including endonuclease activity (GO:0004519), endopeptidase activity (GO:0004175), DNA metabolic processes (GO:0006259), and kinesin complex (GO:0005871), among others (Supplementary Table 3). It should be noted that due to the low representation of genic CpGs in our RRBS data, we are missing what could be key information about CpG methylation in gene regions. Therefore, future studies should further examine CpG methylation across gene regions using a whole genome bisulfite sequencing approach. However, this dataset reveals what could be important genes and functional categories associated with DNA methylation in sea urchins. In summary, we identified multiple genic CpGs that varied by environmental context, particularly maternal, with functions indicative of DNA integration, modification and turnover of proteins and nucleic acids as well as cytoskeletal components.

To further examine potential connections between environmentally induced changes in DNA methylation and plastic phenotype, such as spicule length, we also sequenced transcriptomes of prism stage embryos. RNA sequencing yielded an average of 29,722,500 reads per sample across 12 samples (Supplementary Table 1). After size and quality trimming, an average of 29,713,377 reads per sample remained. Trimmed reads were mapped to the Strongylocentrotus purpuratus genome (V3.1) with an average mapping efficiency of 62.1%. Mapped reads were converted to gene counts across 30,284 gene meta-features with an average of 7,897,000 counts per sample (Supplementary Table 1). The outlier detection package, arrayQualityMetrics, did not detect any sample outliers so all were retained. Genes with low mean counts (less than 3) were discarded from further analysis, leaving 18,015 genes total.

DESeq2 models revealed significant differentially expressed genes for both maternal and developmental conditions. Differences in maternal environment resulted in 1,811 DEGs while developmental environment resulted in 3,765 DEGs (FDR < 0.05), with 645 DEGs overlapping between the two factors (Supplementary Figure 3). PCoAs of gene counts revealed strong differences in gene expression between each of the four treatments where developmental condition separated variation along PCo1 and maternal condition separated variation along PCo2 (Figure 1B). A permutational multivariate ANOVA revealed that 29.9% (p = 0.001) of the variation in gene expression was explained by developmental condition while 19% (p = 0.002) was explained by maternal condition (Figure 1D). The interaction between the two factors explained 6.6% of the variation, although this was not significant. Overall, a large proportion of the prism-stage transcriptome was responsive to different environmental conditions; this is not likely due to developmental delay as staging was performed using developmental cues not body size.

To characterize environmentally responsive transcriptional changes, we performed a weighted gene co-expression network analysis (WGCNA). WGCNA revealed 12 gene expression modules with all showing significance for at least one experimental treatment or physiological trait (Supplementary Figure 5). Each module contained between 2,676 (turquoise) and 84 genes (royalblue) and the majority of modules showed significant enrichment of GO categories. Significant modules were binned as either strongly representative of changes due to developmental condition (turquoise, pink, red, black, and yellow) (Figure 2) or maternal condition (lightyellow, cyan, royalblue, lightgreen, blue, and red) (Figure 3). There was a trend for maternal modules to have a higher percentage of genes with DMCpGs in them, although this trend was not significant (Supplementary Figure 5). The red module represents genes that show significant correlations for both experimental treatments and therefore represents interaction genes, however, there was no significant enrichment of GO terms within this module (Supplementary Table 4).

Developmental modules include those with genes that were up-regulated when embryos were reared in upwelling as compared to non-upwelling conditions (turquoise and pink, Figure 2) and those with genes that were down-regulated when reared in upwelling as compared to non-upwelling (black, Figure 2). The turquoise and pink modules show enrichment of GO categories associated with RNA processing, methylation and metabolic processes, as well as cellular response to stress and growth. In contrast, genes in the black module show larvae reared under upwelling conditions down-regulated multiple neurological and signaling processes, suggesting these physiological processes are suppressed under upwelling stress. Calcium ion binding genes were also reduced when larvae were reared in upwelling conditions, which correlate strongly with the noticeable reduction in spicule length of these embryos (Supplementary Figure 4). The yellow module showed no significant enrichment of GO terms (Supplementary Table 4).

Of the modules associated with differences in maternal environment, only two of them showed significant enrichment of GO terms, the cyan and blue modules (Figure 3). The cyan module represents genes down-regulated in larvae when mothers were reared in upwelling conditions and includes changes in DNA replication and metabolic processes (Figure 3B). The blue module represents genes upregulated in larvae when mothers were reared in upwelling conditions and also includes GO categories associated with DNA metabolism, translation and hydrogen ion transport.

Due to environmental effects on DNA methylation, transcription, and spicule length, we further examined patterns that might indicate causality. Percent methylation values per gene were associated with overall gene expression. There was a strong significant correlation between mean gene expression counts and percent methylation (Adjusted R² = 0.1836, p_lm < 2.2e–16) (Supplementary Figure 7). A total of 3,359 genes had >5 sequenced CpG sites and gene expression counts >3. However, there was little association between differences in GS-CpG methylation and differential expression of genes (Figures 4, 5). Differential gene expression was evident in both genes with low methylation as well as high methylation (Figure 4); however, DEGs showed higher log₂FC differences amongst lowly methylated genes, particularly when considering differences in maternal treatment (Figure 4B). Genes with DMCpGs were almost exclusively those with high methylation, regardless of experimental treatment (Figure 4). There was a significant negative relationship between developmental DEGs and percent methylation; lowly methylated genes were more likely to be upregulated when larvae were reared under upwelling conditions (R² = 0.11, p_lm < 0.001, Figure 4A).

There was some overlap between genes with DMCpGs and DEGs from the DESeq2 models in response to the treatments (Figure 5). Of the 136 genes with maternal DMCpGs, 12 were differentially expressed with respect to maternal environment and 22 were differentially expressed with respect to developmental environment (Figure 5 and Supplementary Table 5). Of the 49 genes with developmental DMCpGs, 5 were differentially expressed with respect to maternal environment and 12 were differentially expressed with respect to developmental environment (Figure 5 and Supplementary Table 5). Overall, very few individual genes exhibited both differential expression and DMCpGs, suggesting it is unlikely that changes in DNA methylation at the level of gene bodies impart regulatory control on transcription, at least across the time frames examined in this study.

The transcriptome is highly sensitive to environmental differences, ultimately changing organismal traits from the molecular level up to higher-level processes such as behavior. Environmental perturbations can induce substantial changes in the transcriptome within short time frames (on the scale of minutes). These fast reaction times are necessary to alter cellular machinery in response to change. DNA methylation, on the other hand, is relatively stable and while able to change in response to the environment, does so at a lesser magnitude than gene expression and across longer time scales. Therefore, in the context of our experimental design, we see the biggest differences in larval CpG methylation when mothers were conditioned in different environments for 4 months, while the developmental treatment was only ∼72 h. Because of the stability and proposed functions of DNA methylation in gene regions, modulating CpG methylation during periods of short-term environmental differences is not likely to be an evolutionarily favored response, since short-term stress might not be indicative of persistent environmental change. Longer-term chronic environmental differences, such as maternal conditioning in the current experiment, may be more predictive of the future environment, so modulation of DNA methylation in response to this may be an adaptive response that occurs across longer time scales and generations.

In this experiment, we chose to use conditions representative of upwelling, which is variable and periodic in the Santa Barbara Channel. This design employed was intentionally simplified from the in situ environmental conditions in an attempt to identify the potential for environmentally induced methylation and gene-expression changes and query any mechanistic relationship between methylation changes driven by the maternal and developmental environment and any corresponding change in transcription. Through field experiments, we know that even variable in situ conditions can drive maternal effects (Hoshijima and Hofmann, 2019). Therefore, going forward it will be important to investigate whether offspring DNA methylation patterns are modulated in response to non-static maternal conditioning or if the induction of methylation changes is dependent upon the predictability of the environmental changes (Burgess and Marshall, 2014).

Mechanisms connecting DNA methylation with transcription are less well understood in invertebrates with sparsely methylated genomes as opposed to vertebrate systems where promoter methylation is able to directly alter transcription and downstream phenotypes (reviewed in Li and Zhang, 2014). Generally, gene body methylation is positively correlated with gene expression levels across invertebrate taxa (Zemach et al., 2010) including, bivalves (Gavery and Roberts, 2013; Wang et al., 2014), ascidians (Suzuki et al., 2007), and social insects (Bonasio et al., 2012; Sarda et al., 2012; Wang et al., 2013; Patalano et al., 2015). In corals specifically, it has been suggested that DNA methylation fine tunes gene expression in response to environmental change via reduction in spurious transcription (Li et al., 2017; Liew et al., 2018), genome-wide balance of gene expression (Dixon et al., 2018) and through shaping codon usage over evolutionary time scales (Dixon et al., 2016). Similar studies have been shown in oysters as well (Gavery and Roberts, 2010; Roberts and Gavery, 2012; Wang et al., 2014; Jiang et al., 2015; Saint-Carlier and Riviere, 2015). While these functions are not mutually exclusive, consensus on how actively DNA methylation regulates transcription in response to environmental change is not well agreed upon across many invertebrate systems. That being said, our data is suggestive that DNA methylation changes in genes rarely correspond to gene expression on a gene by gene level, despite broad genome-wide patterns: similarly enriched GO terms are found between DEGs and genes with DMCpGs associated with different maternal environments, highly methylated genes are more likely to be highly transcribed, and genes with DMCpGs are always highly methylated.

There is contention as to whether gene expression or DNA methylation imparts control on the other. While the overall magnitude of gene expression responses to environmental conditions is larger, there was a stronger influence of maternal conditioning on DNA methylation than of the abiotic conditions under which the embryos developed. This might suggest that maternally mediated changes in DNA methylation may translate to differential expression in offspring. However, we found minimal overlap between maternally driven differences in DNA methylation and transcription of those genes in their offspring; the majority of genes with DMCpGs by maternal condition were not differentially expressed. However, genes with DMCpGs showed functional enrichment of DNA metabolic processes and DNA integration (Supplementary Table 3), similar to terms enriched in modules showing strong correlation with maternal conditions (Figure 3 and Supplementary Table 4). This is evidence that while there are environmentally inducible changes in both CpG methylation and gene expression, they appear to be only broadly linked. This suggests that environmentally induced changes in DNA methylation may play a different functional role on a cellular level than regulating transcription directly as is observed in vertebrate systems, potentially through controlling spurious transcription. However, it should be noted that the reduced representation approach limits our ability to examine methylation patterns in all genes and a full genomic picture connecting transcription and DNA methylation is warranted.

While the focus of our study was to examine the role of DNA methylation in response to environmental conditions, it is worth noting that DNA methylation is also believed to be necessary to ensure successful development in invertebrates (Regev and Lamb, 1998; Riviere et al., 2013). Indirect evidence of this in sea urchins includes observations that DNA methylation of broad functional groups and developmentally regulated genes vary by life history stage in S. purpuratus (Fronk et al., 1992; Strader et al., 2019). Additionally, DNA methyltransferase activity has been shown to change dramatically during the early development of the sea urchin Sphaerechinus granularis (Tosi et al., 1995). Treatment with 5-azacytidine, a nuclear methylation inhibitor, at any stage prior to blastula inhibits development beyond the blastula stage in sea urchin embryos, Paracentrotus lividus and Sphaerechinus granularis, causing embryos to arrest at the blastula stage though some continue to develop spicules with no other morphological change (Crkvenjakov et al., 1970; Maharajan et al., 1986). Although details regarding the specific role that DNA methylation plays in sea urchin development are generally lacking, there is evidence that it is critical to the development of other invertebrates such as the oyster, Crassostrea gigas (Riviere et al., 2013), the wasp, Nasonia vitripennis (Zwier et al., 2012), and the ascidian, Phallusia manzmilata (Maharajan et al., 1986). Thus, in order for development to proceed successfully, we might expect certain limitations on how much DNA methylation can change in developing embryos and larvae. Our previous work showed minimal differential methylation in S. purpuratus as a function of developmental conditions (Strader et al., 2019), although the experiment included developmental environments that only varied by pCO₂ level. In this study, the combined temperature and pCO₂ levels of our non-upwelling versus upwelling developmental treatments resulted in differential methylation, but to a lesser extent than the effects of maternal environment on differential methylation. Therefore, some alterations in DNA methylation can occur in response to external factors in spite of potential constraints on methylation during early development, although this likely varies by the multitude and intensity of the environmental stressor(s).

The blue module shows an upregulation of genes associated with hydrogen ion transmembrane transporters and hydrogen transport when larvae had mothers that were reared in upwelling conditions. This suggests that there is a strong maternal effect of these genes, potentially priming these larvae to express more H⁺ transporters in the high pCO₂ environment associated with the upwelling treatment. While extracellular tissue in larval urchins maintains the pH of the environment, there is internal control of pH within intracellular compartments, including those surrounding the principal calcifying cells. This internal control is mediated by transporters that regulate the concentrations of H⁺ and HCO3- ions on either side of the membrane (Evans et al., 2013). Several studies on gene expression responses to high pCO₂ in larval urchins and other marine invertebrates find either no differential expression or downregulation of hydrogen ion transporters (Todgham and Hofmann, 2009; Moya et al., 2012; Evans and Watson-Wynn, 2014). These studies, however, were all performed by exposing early life-history stages to different pCO₂ environments, similar to the developmental treatment in the present study, which found no differential expression of H⁺ transporters in response to upwelling developmental conditions (Figure 2 and Supplementary Table 5). This suggests that modulation of H⁺ transporters is not likely a fast response to environmental change, as in, within the course of embryonic and larval development. However, longer-term exposure of adult urchins to high pCO₂ and low temperatures (upwelling) can impart an upregulation of these genes that are important in maintaining the balance between H⁺ and HCO3- ions in offspring. The upregulation of these genes, however, does not seem to generate a benefit to the larvae when exposed to upwelling conditions themselves with regards to larval spicule length (Supplementary Figure 4) and expression of cellular stress response genes (Figure 2).

Black and cyan modules showed enrichment of genes associated with DNA metabolic processes, macromolecule biosynthetic processes and DNA integration, suggesting overall differential regulation of these processes when mothers are exposed to different conditions. Functional enrichment of genes with DMCpGs, the majority differentially methylated with respect to maternal environment, also reveals enrichment of DNA integration, macromolecule biosynthetic processes, and DNA metabolic processes, among others (Supplementary Table 3). Despite similarities in GO terms, maternal modules did not contain more genes with DMCpGs than developmental modules (Supplementary Figure 5) and there were very few genes that exhibited both DMCpGs and DEGs by maternal condition (12 genes, Supplementary Table 5). The result that maternal responses to upwelling imparted differential methylation and expression in genes with similar functionality suggests that modulating the activity of these genes may be a critical long-term response to upwelling conditions. While there are caveats in extrapolating GO categories from only a subset of genes with sequenced CpGs, in Aiptasia sea anemones, similar GO categories are enriched in DM genes, including “DNA-dependent DNA replication” and “DNA integration” (Li et al., 2017), the same top GO category enriched for genes with DMCpGs in this dataset (Supplementary Table 3). DNA integration is the process by which one segment of DNA is inserted into another, one example being transposable elements. Epigenetic and transposable elements strongly interact with each other and both can lead to changes in phenotypes and genotypes in response to stress (Rey et al., 2016). For example, epigenetic components are key in repressing transposable element activity. Therefore, it is possible that interactions between DNA methylation and transposable elements can drive rapid changes in phenotypes in response to global change across short time scales in taxa other than vertebrates and plants (Rey et al., 2016). Our results, as well as others reporting differential methylation of genes associated with DNA integration, posit that more research into potential effects of transposable elements on rapid phenotypic change is warranted in marine systems exposed to a global change multi-stressor scenario. In addition, DNA methylation can impact phenotypes by changing the mutability of gene sequences, which has been suggested to impact codon bias evolution (Dixon et al., 2016). These mechanisms might be better representative of down-stream effects of environmentally induced changes in DNA methylation than directly changing transcription in invertebrates with sparsely methylated genomes.

The developmental environment had a more marked impact on spicule length and the transcriptome than did the adult environment during gametogenesis (Figure 1 and Supplementary Figure 4). This is in contrast to a previous study in S. purpuratus in which the maternal environmental conditions appeared to have a greater impact on the offspring phenotype and transcriptome than the conditions during development (Wong et al., 2018). This dissimilarity could be in part due to differences in developmental stage, as organismal responses to the environment have been shown to vary greatly by developmental stage (Kurihara, 2008; Ross et al., 2011; Byrne, 2012). Wong et al. (2018) examined the body size and transcriptomic patterns of the gastrula stage, during which major developmental landmarks distinguish this from later larval stages. Another noteworthy distinction is that in Wong et al. (2018), the embryos were raised under two different conditions that only varied by a single abiotic factor: pCO₂ level. In the study presented here, combinations of both temperature and pCO₂, intended to reflect ecologically relevant upwelling or non-upwelling conditions, were manipulated during early development so as to generate similar or reciprocal treatments as the adults experienced. Thus, the dominance of phenotypic and transcriptomic plasticity as a result of developmental conditions may be a result of the inclusion of temperature as a factor during development, or due to the multi-stressor nature of including combinations of different temperature and pCO₂ levels.

Offspring that were reared in upwelling conditions exhibited phenotypic and transcriptomic patterns associated with a response to stress (Figure 2 and Supplementary Table 4). In concert with this stress response, there was an up-regulation of numerous RNA modification processes as well as processes involved in translation and potentially post-translational modifications (protein polymerization). For example, one gene in the turquoise module, also with a significantly differentially methylated CpG (Figure 5B), Sp.Dusp6.7.9, is known to be involved with protein dephosphorylation. Finally, there was also upregulation of genes involved in carbohydrate and lipid transport. These results strongly suggest larvae reared in upwelling conditions exhibit a stress response and a potential reallocation of energy resources, mostly likely to maintain calcification and growth. Larvae reared in upwelling conditions also exhibit a significant reduction in spicule length compared to those reared in non-upwelling conditions (Supplementary Figure 4). Differences in RNA metabolic and processing pathways suggest an attempt to modulate cellular functions in response to upwelling stress. In addition, the majority of DMCpGs were hypermethylated in the upwelling condition compared to the non-upwelling condition, which is a similar response to other organisms exposed to pH stress (Liew et al., 2018). Overall, these results show clear compensatory mechanisms in response to upwelling stress.

Stress as a result of developing under the simulated upwelling treatment could be due to low pH, low temperature, or the combination of both. In regards to low pH environments, sea urchin species have generally been shown to display reduced skeletal growth and calcification (Kurihara and Shirayama, 2004; Clark et al., 2009; Dupont et al., 2010; Byrne and Przeslawski, 2013), similar to what we found in this study. The increased expression of genes related to calcium homeostasis has been suggested as a means by which developing sea urchin larvae are able to cope and maintain skeletogenesis in low pH environments (Evans et al., 2013; Evans and Watson-Wynn, 2014). Thus, failure to upregulate calcium-related genes could result in decreased skeletogenesis. Indeed, larvae in this study exhibited a relative downregulation of calcium ion binding genes, a similar pattern to that observed by Todgham and Hofmann (2009) in S. purpuratus larvae raised under elevated pCO₂ levels. In regards to temperature effects, sea urchin larvae across a variety of species and latitudes generally exhibit greater growth and calcification when raised under warmer temperatures, which occasionally provides a compensatory effect that helps offset the negative effect of low pH (Sheppard Brennand et al., 2010; Byrne, 2012; Byrne and Przeslawski, 2013; Wangensteen et al., 2013). Throughout upwelling events, however, no such compensatory effect is expected because, in addition to higher pCO₂, cold temperatures are associated with upwelling seawater. Therefore, the down-regulation of calcium ion binding genes and the significant reduction of spicule length under upwelling conditions align with our expectations for larvae developing under low pH and/or low temperature conditions.

Epigenetic inheritance, in which epigenetic marks are intergenerationally transmitted through the germline, is a major study question within the field of environmental epigenetics (Verhoeven et al., 2016; Eirin-Lopez and Putnam, 2019), as it may contribute to rapid acclimation and resilience to changing environments, thereby influencing the fitness landscape and resulting evolutionary processes (Bonduriansky and Day, 2009; Bonduriansky et al., 2012). Although this study does not provide direct evidence of this, the pattern in which maternal conditions had a greater influence on differential methylation between larvae than developmental conditions supports the potential for epigenetic inheritance. In this case, epigenetic marks in the form of DMCpGs may have been incorporated by the mothers into their eggs, which then persisted in their offspring during early development despite differences in the environmental conditions under which the embryos and larvae were raised. While further study is required to determine if these epigenetic changes are truly integrated into the germline and are capable of remaining across generations, this study provides support for the potential of epigenetic inheritance in this system.