Gene Expression Patterns of Red Sea Urchins (Mesocentrotus Franciscanus) Exposed to Different Combinations of Temperature and pCO2 During Early Development

doi:10.21203/rs.3.rs-36920/v1

Download PDF

Research article

Gene Expression Patterns of Red Sea Urchins (Mesocentrotus Franciscanus) Exposed to Different Combinations of Temperature and pCO2 During Early Development

https://doi.org/10.21203/rs.3.rs-36920/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 07 Jan, 2021

Read the published version in BMC Genomics →

You are reading this older preprint version

Read the latest preprint version →

Background: The red sea urchin Mesocentrotus franciscanus is an ecologically important kelp forest herbivore and an economically valuable wild fishery species. To examine of how M. franciscanus responds to its environment on a molecular level, differences in gene expression patterns were observed in embryos raised under combinations of two temperatures (13 °C and 17 °C) and two pCO₂ levels (475 matm and 1050 matm). The transcriptomic responses of the embryos were assessed at two developmental stages (gastrula and prism) in light of previously described plasticity in body size and thermotolerance under these temperature and pCO₂ treatments.

Results: Although transcriptomic patterns primarily varied by developmental stage, there were pronounced differences in gene expression as a result of the treatment conditions. Temperature and pCO₂ treatments led to the differential expression of genes related to the cellular stress response, transmembrane transport, metabolic processes, and the regulation of gene expression. Temperature had a greater influence on gene expression than pCO₂, and may have contributed to positive effects of temperature on body size and thermotolerance at the prism stage. On the other hand, a relatively muted transcriptomic response to pCO₂ may have permitted the stunting effect of elevated pCO₂ on embryo body size.

Conclusions: M. franciscanus exhibited both transcriptomic and phenotypic plasticity in response to temperature and pCO₂ stress during early development. As climate change continues, red sea urchins may benefit from moderate ocean warming, whereas they will be negatively affected by ocean acidification. Present-day pCO₂ conditions that occur due to coastal upwelling may already be detrimental to populations of M. franciscanus.

Epigenetics & Genomics

red sea urchin

Mesocentrotus franciscanus

RNA-seq

transcriptomics

early development

climate change

warming

ocean acidification

The red sea urchin Mesocentrotus franciscanus (A. Agassiz, 1863) is an ecologically and economically valuable species found along the Pacific Coast of western North America [1]. In subtidal areas, especially within kelp forests, these echinoderms are herbivorous ecosystem engineers that can shape the flow of resources within marine habitats [2]. Overgrazing by M. franciscanus, often in combination with overgrazing by the purple sea urchin Strongylocentrotus purpuratus, can lead to the formation of urchin barrens in which macroalgal communities are severely reduced or depleted [3, 4]. Red sea urchins also function as prey to animals at higher trophic levels, including spiny lobsters and sea otters [5–7]. In addition to its removal by natural predators, M. franciscanus is widely collected as a lucrative wild fishery species. Fisheries in Mexico, the United States, and Canada harvest M. franciscanus for their gonads (i.e., roe) that supply domestic markets as well as international exports, principally to Japan [8, 9]. Over the last five years (2015–2019), the annual revenue reported from M. franciscanus fisheries across the states of California, Oregon, and Washington averaged over $7.1 million USD/year, far more than all other echinoderm fishery species combined [10].

Given the considerable ecological and economic importance of M. franciscanus, determining how this species will be affected by continuing environmental change in coastal oceans remains an overlooked and critical area of research [11]. Due to their habitat and life history, these urchins are threatened by climate change impacts [12] such as ocean warming, which may include sudden and extreme marine heat waves [13, 14], and ocean acidification, which may amplify the low pH conditions that episodically occur in upwelling regions [15]. Furthermore, these urchins may be particularly vulnerable to stress during early development. Although planktonic embryological and larval stages of echinoids are capable of exhibiting vertical migration [16, 17], they are likely less capable of finding refuge from stressful conditions than their benthic adult counterparts. There is also evidence that many organisms are most vulnerable to environmental stress early in their life history [18–21]. Both lethal and sublethal effects that occur during early development or that carry over into later life stages will negatively affect the recruitment necessary to support future populations [22, 23]. A limited number of studies have investigated how M. franciscanus responds to temperature or pCO₂ stress [24–26], and even fewer have done so within a multi-stressor context [27]. Here, both temperature and pCO₂ conditions were manipulated in a laboratory setting to investigate the gene expression patterns of M. franciscanus during its early development. To the best of our knowledge, this is the first study to use RNA sequencing (RNA-seq) to examine the M. franciscanus stress response.

In this study, M. franciscanus embryos were raised under a combination of two temperatures (13 °C or 17 °C) and two pCO₂ levels (475 µatm or 1050 µatm) that reflect current and future ocean conditions in their natural habitat [15, 28–30]. This generated four different treatment combinations: 1) 17 °C and 1050 µatm pCO₂, 2) 17 °C and 475 µatm pCO₂, 3) 13 °C and 1050 µatm pCO₂, and 4) 13 °C and 475 µatm pCO₂. Gene expression patterns were assessed at both the gastrula and prism embryo stages. Physiological assessments from this experiment, including body size and thermotolerance, have been previously reported [31]. Here, we describe the effects of the temperature and pCO₂ treatments at the molecular level and whether they relate to observations made at the level of the organism. Temperature elicited a robust transcriptomic response at both developmental stages. Gene expression analyses indicated that the warmer temperature (i.e., 17 °C) induced cellular stress, amongst other processes. Nonetheless, the transcriptomic response may have contributed to no effect of temperature on body size at the gastrula stage, increased growth under warmer temperatures at the prism stage, and slightly greater thermotolerance. In contrast, a comparatively muted transcriptomic response to the pCO₂ treatment may be responsible for smaller-sized embryos under elevated pCO₂ levels (i.e., 1050 µatm). Overall, we examined a valuable fishery species that is capable of dramatically shaping coastal ecosystems, and determined that during early development M. franciscanus exhibits both transcriptomic and phenotypic plasticity in response to two climate change-related stressors. This study provides much needed insight into a species that is important for many local fisheries on the Pacific coast of North America, facilitating our understanding of how M. franciscanus development is affected by current ocean conditions, as well as our predictive capacity of how this species will respond to future ocean change.

Summary statistics and overview of RNA-seq

The samples used for RNA-seq were generated from triplicate cultures of embryos raised at each of the four combined temperature and pCO₂ treatments (i.e., 12 total cultures). Each sample was collected as a pool of 5000 embryos from each of the 12 cultures at both the gastrula and prism stages during development to produce a total of 24 samples used for RNA extractions and library preparation. Sequencing of the 24 libraries yielded a total of 728,782,735 100-bp single reads. After quality trimming, an average of 30.3 ± 1.3 million reads per library remained. FASTQC reports [32] of trimmed sequences showed high sequence quality (> 30) with limited adapter contamination or presence of overrepresented sequences. Per-library mapping efficiency to the developmental transcriptome [33] using RSEM [34] was at an average of 52.6%. After filtering the data to those with more than 0.5 counts per million mapped reads across at least three samples, 37,577 sequences remained.

Developmental stage influenced transcriptomic patterns

A principal component analysis (PCA) of sample-to-sample distances showed that differences in gene expression profiles were primarily between the two developmental stages, gastrula and prism (Fig. 1a). Principal Component 1 (PC1) captured the majority of the variance (58.6%) and revealed a clear separation between gastrula and prism stage embryos, while Principal Component 2 (PC2) only captured 4.6% of the variance. Indeed, a permutational multivariate ANOVA across all 24 samples with developmental stage, temperature treatment, and pCO₂ treatment as fixed factors, revealed that developmental stage explained 57.8% of the variance (p < 0.001) (Fig. 1b). Because we have previously explored the differences in gene expression patterns across M. franciscanus during early development [33] and it is not the main focus of the current study, from here onward we report separate gene expression analyses for the gastrula and prism stages.

Temperature and pCO₂ treatments affected gastrula gene expression

Separate PCA plots were generated for the gastrula and prism stages. At the gastrula stage, we generally observed that both temperature and pCO₂ treatments appeared to drive differences in gene expression patterns across samples. A PCA of only the gastrula stage showed that replicate samples grouped together (Fig. 1c). Here, PC1 captured 23.1% of the variance and appeared to separate samples by temperature treatment (17 °C versus 13 °C). PC2 captured 11.2% of the variance and appeared to separate the samples by pCO₂ treatment (1050 versus 475 µatm).

A permutational multivariate ANOVA revealed that at the gastrula stage, 19.2% of the variance was explained by temperature treatment (p < 0.001) (Fig. 1d). Differential expression (DE) analyses conducted in limma [35] identified differentially expressed genes (relatively up- and down-regulated) between gastrula raised under different temperature treatments. A total of 4636 genes were significantly up-regulated in embryos raised at 17 °C relative to embryos raised at 13 °C (p < 0.05) (Fig. 2a). Gene ontology (GO) analyses were performed following the GO_MWU package [36] in R and identified significantly enriched terms (p > 0.05) across molecular function (MF), biological process (BP), and cellular component (CC) GO categories for genes differentially expressed due to temperature (Fig. 3a, Additional file 1). GO categories significantly enriched among up-regulated genes included DNA recombination, DNA metabolic process, unfolded protein binding, and G protein-coupled receptor (Fig. 3a). A total of 1735 genes were down-regulated in embryos raised at 17 °C relative to embryos raised at 13 °C (Fig. 2a). GO categories enriched with down-regulated genes included regulation of gene expression, chromatin organization, negative regulation of biological process, and ion binding (Fig. 3a).

The pCO₂ treatment also affected gene expression patterns at the gastrula stage, explaining 12.5% of the observed variance (p = 0.034) (Fig. 1d). GO analyses identified significantly enriched terms (p > 0.05) for genes differentially expressed due to pCO₂ (Fig. 3b, Additional file 1). Only 49 genes were up-regulated when comparing the 1050 µatm to the 475 µatm pCO₂ treatment at the gastrula stage (Fig. 2a). GO categories enriched with up-regulated genes included macromolecule catabolic process, ion binding, active transmembrane transporter, ATPase, and heat shock protein binding (Fig. 3b). A total of 202 genes were down-regulated in embryos raised at 1050 µatm relative to embryos raised at 475 µatm (Fig. 2a). Enriched GO categories included biosynthetic processes, nucleic acid binding, and ncRNA metabolic process (Fig. 3b). The interaction between temperature and pCO₂ factors explained 7.4% of the variance observed at the gastrula stage, but the interaction was not significant (p = 0.419) (Fig. 1d).

Temperature treatment affected prism gene expression

Similar to the gastrula stage, the PCA of only the prism stage showed a separation of samples by experiment treatment with sample replicates grouping together (Fig. 1e). Along the PC1 axis, which captured 23.1% of the variance, there was a clear separation of samples by treatment temperature (17 °C versus 13 °C). Along the PC2 axis, which captured 10.2% of the variance, there appeared to be a separation of samples by pCO₂ treatment level (1050 versus 475 µatm). A permutational multivariate ANOVA revealed that at the prism stage, 22.8% of the variance was explained by the temperature treatment (p < 0.001) (Fig. 1f). DE analysis showed a total of 4132 genes were up-regulated in embryos raised at 17 °C relative to those raised at 13 °C (Fig. 2b). GO analysis revealed GO categories enriched in these up-regulated genes, which included oxidoreductase, response to oxidative stress, ion binding, ion transmembrane transporter (Fig. 4, Additional file 1). A total of 4286 genes were down-regulated in embryos raised at 17 °C relative to those raised at 13 °C (Fig. 2b). Enriched GO categories included regulation of gene expression, RNA methyltransferase, histone modification, cellular response to DNA damage stimulus, and response to unfolded protein (Fig. 4). The pCO₂ treatment only explained 9.6% of the variance at the prism stage and was not significant (p = 0.084) (Fig. 1f). In fact, only 60 genes were down-regulated and 3 genes were up-regulated when comparing the 1050 µatm to the 475 µatm pCO₂ treatment at the prism stage (Fig. 2b). Additionally, the interaction between temperature and pCO₂ factors explained only 7.7% of the variance and was not significant (p = 0.296) (Fig. 1f).

Gene expression modules were correlated with treatments and physiology

A Weighted Gene Co-Expression Network Analysis (WGCNA) [37] identified genes that were similarly expressed and determined their correlation to each trait of interest (i.e., gastrula temperature treatments, gastrula pCO₂ treatments, gastrula body size, prism temperature treatments, prism pCO₂ treatments, prism body size, and prism thermotolerance). The filtered, normalized and voom-transformed data consisting of 37,577 genes were assigned into ten module eigengenes, while 32 genes remained unclustered and unassigned (represented by the grey module) (Fig. 5). Each module eigengene was related to the traits of interest to generate eigengene networks with positive or negative correlation values ranging from 1 to -1. All ten module eigengenes were significantly correlated to multiple traits of interest (Fig. 5).

WGCNA modules significantly correlated with gastrula traits related to experiment conditions (p < 0.05) showed that five modules were shared across temperature and pCO₂ treatments (pink, red, turquoise, black, and green), while three modules were unique to the temperature treatment (blue, magenta, and brown) and two modules were unique to the pCO₂ treatment (greenyellow and yellow) (Fig. 5). Gene expression modules correlated to gastrula body size, as measured by length in millimeters (mm), were most similar to modules correlated to the lower gastrula pCO₂ treatment (475 µatm). Specifically, gastrula length and 475 µatm treatment traits shared correlation patterns with six modules (greenyellow, pink, red, turquoise, yellow, and black) (Fig. 5).

In contrast, modules correlated to prism traits were largely separated by experimental treatments (Fig. 5). Prism temperature and pCO₂ treatments only shared significant correlations with two modules (pink and red), while six modules were unique to the temperature treatment (blue, greenyellow, magenta, turquoise, black and green) and two modules were unique to the pCO₂ treatment (brown and yellow). The prism body size trait exhibited nearly the same module correlation patterns as the higher prism temperature treatment of 17 °C, sharing the same positive and negative correlations to the same seven modules (blue, greenyellow, magenta, pink, turquoise, black, and green). The thermotolerance trait was even more similar to the 17 °C treatment, sharing the same positive and negative correlations to all eight modules, including the red module.

GO analyses identified significantly enriched terms (p > 0.05) across MF, BP, and CC categories for only three of the modules: blue, turquoise, and black (Fig. 6, Additional file 2). The blue module, composed of 16,164 genes, was positively correlated to the higher temperature treatment (17 °C), and negatively correlated to the lower temperature treatment (13 °C), for both the gastrula and prism stages (Fig. 5). The blue module was also positively correlated to gastrula body size, prism body size, and prism thermotolerance (Fig. 5). GO analysis of the blue module provided GO terms related to response to oxidative stress, active transmembrane transporter, ion transport, and ATP metabolic process (Fig. 6). The turquoise module, composed of 17,316 genes, was positively correlated to the 17 °C gastrula treatment, the 475 µatm pCO₂ gastrula treatment, gastrula body size, and the 13 °C prism treatment (Fig. 5). Likewise, it was negatively correlated with the 13 °C gastrula treatment, the 1050 µatm pCO₂ gastrula treatment, the 17 °C prism treatment, prism body size, and prism thermotolerance. GO terms within the turquoise module were related to regulation of gene expression, histone binding, chromatin organization, RNA modification, and cellular response to DNA damage stimulus (Fig. 6). The black module, which contained 372 genes, was positively correlated with the 13 °C gastrula treatment, the 1050 µatm pCO₂ gastrula treatment, the 17 °C prism treatment, prism body size, and prism thermotolerance (Fig. 5). Conversely, it was negatively correlated with the 17 °C gastrula treatment, the 475 µatm pCO₂ gastrula treatment, gastrula body size, and the 13 °C prism treatment. GO analysis revealed GO terms within the black module included those related to ATPase, ion binding, and cilium movement (Fig. 6).

In this study, we examined how the gene expression patterns of gastrula and prism embryos varied by the developmental temperature and pCO₂ conditions under which they were raised. We also assessed whether the transcriptomic results aligned with the morphometric and physiological results previously reported in [31]. Although both temperature and pCO₂ can influence rates of sea urchin development [38, 39], any potential differences in developmental timing should not have impacted the results of this study because samples were collected based on developmental progression to the desired embryonic stages as detailed in the Methods, rather than by hours post-fertilization. Overall, we found that while transcriptomic patterns varied by developmental stage, temperature had a dominant effect on changes in gene expression while pCO₂ elicited a more subtle transcriptomic response that was largely limited to the gastrula stage. Experimental conditions impacted genes related to the cellular stress response, transmembrane transport, metabolic processes, and the regulation of gene expression.

In terms of experimental design, embryos were obtained by evenly pooling eggs from five females and fertilizing them with sperm from a single male to produce all full or half siblings. Admittedly, there are caveats to this approach. The results presented here may only be representative of a small subset of the population, or it may be driven by the quality of the particular male selected to fertilize the eggs. For instance, in the purple sea urchin S. purpuratus, genetic variation has been shown to influence transcriptomic responses to temperature and pCO₂ stress during early development [40, 41]. We therefore recommend that additional studies be performed within other M. franciscanus populations and with multiple male-female crosses to determine if our results are unique to this study. Nevertheless, this approach was implemented in an effort to limit genetic variability and male-female interactions that may have otherwise confounded the molecular results of this study.

All samples used for RNA extractions were each composed of a pool of 5000 individuals and should thus represent the same mixture of genotypes. Therefore, we do not expect differences in gene expression patterns to be due to genetic variability between embryo cultures, particularly because a low incidence of mortality was observed during the experiment, although it was not directly measured in this study. In the absence of selection, the observed variability in gene expression, body size, and thermotolerance between embryos raised under different experimental treatments reflect plasticity exhibited by M. franciscanus during its early development. We discuss this plasticity, and how it may relate to embryo performance under different conditions that M. franciscanus are likely to experience in their natural environments currently and in the future under ocean change scenarios.

Transcriptomic patterns varied by developmental stage: General patterns

Developmental stage (gastrula or prism embryos) was the primary factor driving differences in gene expression patterns across samples (Fig. 1a and 1b). In a past study, we raised cultures of M. franciscanus embryos in a single laboratory environment that mimicked average, non-stressful conditions in situ (i.e., 15 °C and 425 µatm pCO₂) and documented significant transcriptomic differences between gastrula and prism stages [33]. Therefore, there are many alterations in gene expression between these stages that occur as a result of development and are independent of differences in environmental temperature and/or pCO₂ conditions. This is also evident in Fig. 1a in which gastrula samples do not cluster with prism samples that share the same experimental treatment.

Because comparing gastrula versus prism gene expression patterns was not a goal of this study, no direct differential expression analyses were performed between stages. Nevertheless, embryos at each developmental stage exhibited different transcriptomic responses to temperature and pCO₂ treatments. For instance, at the gastrula stage, many more genes were up-regulated than down-regulated in 17 °C relative to 13 °C, whereas at the prism stage, a similar number of genes were relatively up- and down-regulated between the two temperature treatments (Fig. 2). Additionally, the pCO₂ treatment explained a significant amount of variance in gene expression in gastrula embryos, but not later at the prism stage (Fig. 1d and 1f). The difference in transcriptomic response by developmental stage was also evident in the WGCNA, in which the gastrula and prism temperature and pCO₂ treatments were significantly correlated to many different module eigengenes (Fig. 5). Similarly, the morphometric response to temperature and pCO₂ treatments varied by stage, in which only pCO₂ affected gastrula by reducing body size under elevated pCO₂ conditions (i.e., 1050 µatm) [31]. On the other hand, temperature was the dominant factor at the prism stage, with warmer conditions (17 °C) increasing body size, offsetting the stunting effect of high pCO₂ [31]. The observed patterns between gene expression and body size will be described in greater detail later in the Discussion.

Different life stages are predicted to have different sensitivities to stress [18]. The variability between gastrula and prism stress responses may be explained by a difference in stage-specific vulnerability. During the gastrula stage, the archenteron is formed from invagination of the embryo’s vegetal plate [42], a fundamental process known as gastrulation that is essential for successful development in metazoans [43]. At the prism stage, the embryo develops its digestive tract and skeletal rods, which are vital structures required for the embryos to eventually become feeding, planktotrophic larvae [44, 45]. Accordingly, differences in responses to environmental conditions between these two stages are likely reflective of the distinct processes undergone by these embryos to ensure their continued developmental progression.

The variability between stages could also be due to the timing and duration of exposure to stress. The effects of a stressor can become increasingly deleterious as the length of exposure continues, and organisms not permitted adequate time to recover may exhibit increasingly poor performance. Furthermore, during development there may be negative carry-over effects that persist into later life stages [46, 47]. Alternatively, organisms may acclimate to stressful conditions over time, and are therefore less adversely affected by a stressor following the initial exposure. For example, in the coral Acropora hyacinthus, the immediate transcriptomic response to heat stress was much higher than the transcriptomic response following 20 hours of exposure to warmed conditions [48]. Thus, it remains important to acknowledge that organisms may responds differently to various environmental stressors depending on their life history as well as the timing and duration of the exposure.

Temperature influenced gastrula embryos on a molecular level

Temperature was the dominant factor influencing changes in gene expression at the gastrula stage, explaining 19.2% of the observed variance (p < 0.001) with 4626 genes up-regulated and 1735 genes down-regulated in embryos raised under 17 °C relative to those raised under 13 °C (Figs. 1 and 2). In general, the observed temperature effects on gene expression at the gastrula stage were approximately akin to those reported for the purple sea urchin S. purpuratus [41], whose biogeographical distribution overlaps with that of M. franciscanus. Here, DE analysis and WGCNA revealed that gastrula raised in the higher temperature treatment (i.e., 17 °C) expressed genes associated the cellular stress response. Gastrula embryos of S. purpuratus that were raised under an 18 °C temperature treatment exhibited a comparable cellular stress response by up-regulating genes associated with cellular responses to reactive oxygen species and unfolded proteins [41]. We also found that M. franciscanus gastrula embryos raised under the warmer treatment exhibited transcriptomic patterns indicative of increased transmembrane transport and metabolism, while embryos from the colder treatment appeared to decrease metabolic processes. Similarly, S. purpuratus increased expression of ion channel, cell-cell signaling, and metabolism genes at warmer temperatures [41]. Lastly, temperature appeared to impact the regulation of gene expression in M. franciscnus gastrula embryos, including genes related to epigenetic and epitranscriptomic mechanisms. Temperature also appeared to influence how gene expression was regulated in S. purpuratus gastrula embryos with higher temperatures leading to a down-regulation of genes related to transcription and RNA processing [41].

Despite similarities in how temperature influenced gastrula gene expression, unlike in M. francicanus, there was an effect of temperature on S. purpuratus morphology. Specifically, S. purpuratus gastrula embryos raised at warmer temperatures were significantly smaller in size [41], whereas M. franciscanus gastrula embryos did not significantly differ in size as a result of temperature [31]. This could reflect differences in experimental design between the studies (e.g., treatment temperatures, breeding designs, and urchin collection sites), or it could reflect differences that exist at the species level. Unlike S. purpuratus, the temperature response of M. franciscanus gastrula embryos that occurred at the molecular level was not reflected at the organismal level. We postulate that the transcriptomic differences between gastrula raised at 17 °C and 13 °C served to compensate for direct temperature effects and allowed the embryos to maintain the same size despite the temperature treatments. For the remainder of this section, we explore with greater detail how temperature affected the expression patterns of genes associated with the cellular stress response, transmembrane transport, metabolism, and gene expression regulation in M. franciscanus gastrula embryos.

Cellular stress response

Environmental stress can lead to the production of reactive oxygen species (ROS), which can cause oxidative stress if ROS production exceeds the organism’s antioxidant or damage repair capacity [49, 50]. Oxidative stress, and the response to resulting cellular damage, due to elevated temperatures have been documented across a wide variety of taxa, including algae [51], plants [52], mollusks [53], and fishes [54, 55]. Gastrula embryos raised at 17 °C were positively correlated to the blue WGCNA module that included genes associated with oxidoreductase activity (Figs. 5 and 6). Oxidoreductase enzymes catalyze the transfer of electrons from a reductant to an oxidant [56] and can contribute to the production of ROS [57–59]. The blue module also included genes related to response to oxidative stress and response to ROS (Fig. 6). Thus, at 17 °C, gastrula embryos appear to be responding to oxidative stress at a molecular level.

ROS can cause cellular damage to lipids, proteins, and nucleic acids [59, 60]. Genes related to unfolded protein binding and glycosylation (Fig. 3a) were up-regulated in gastrula embryos raised at 17 °C, indicating a response to protein alteration and damage. The unfolded protein response acts to restore protein folding in the endoplasmic reticulum and reestablish protein homeostasis [61, 62], while glycosylation is important for regulating the structure, function and stability of proteins [63, 64]. DE analysis and WGCNA also provided evidence of DNA damage and repair. This is similar to observations in Acropora corals in which heat stress caused an up-regulation of DNA replication and repair genes [48, 65]. Here, the turquoise WGCNA module, which was positively correlated to the gastrula 17 °C treatment (Fig. 5), included GO terms related to cellular response to DNA damage stimulus as well as DNA replication and repair, such as DNA polymerase, DNA metabolic process, DNA binding, and helicase (Fig. 6). Genes related to DNA metabolic process and DNA recombination were also up-regulated under warmer temperature conditions (Fig. 3a). DNA metabolic processes can include both DNA synthesis and degradation for the purposes of replication and repair. Furthermore, DNA recombination in somatic cells has been identified as a critical mechanism for DNA damage repair [66, 67]. Taken together, 17 °C temperature conditions appear to induce oxidative stress within the gastrula embryos, which undergo stress response mechanisms to combat cellular damage.

Transmembrane transport and metabolism

Gastrula embryos raised under warmer temperatures also increased expression of genes related to transmembrane transport, potentially both within and between cells (i.e., cell-cell communication). For instance, genes related to G protein-coupled receptor, cation channel, cell surface receptor signaling pathway, and vesicle-mediated transport were up-regulated in gastrula embryos raised under 17 °C relative to those raised under 13 °C (Fig. 3a). Additionally, GO terms within the blue WGCNA module, which was positive correlated to the gastrula 17 °C treatment (Fig. 5), included sodium and potassium ion transmembrane transporters, ATP hydrolysis coupled cation transmembrane transport, proton transmembrane transport, and neurotransmitter transport (Fig. 6). Increased transport of materials, particularly ions, across cell membranes may indicate osmoregulation and maintenance of homeostasis. This aligns with reports in juvenile sea urchins of the species Loxechinus albus, in which gene expression alterations under elevated temperatures provided evidence of increased active transmembrane transport of sodium and potassium ions [68].

Active transport of ions across cell membranes is an energetically expensive process that can incur substantial metabolic costs during sea urchin early development [69]. The blue module also contains genes related to lipid metabolic process, peptide metabolic process, carbohydrate metabolic process, and ATP metabolic process (Fig. 6). This up-regulation of metabolic genes is similar to observations in S. purpuratus gastrula embryos raised under warmer temperatures [41]. In contrast, gastrula raised under the lower temperature treatment expressed genes associated with metabolic depression. Specifically, genes up-regulated in gastrula embryos raised under 13 °C relative to those raised under 17 °C included those related to negative regulation of biological process and negative regulation of metabolic process. In this study, metabolic rates of embryos raised under different treatments were not measured at the gastrula stage, but we may expect that, given the effect of temperature on biochemical reaction kinetics, metabolic rate should increase predictably with temperature [70]. Generally, higher metabolic rates have been recorded at warmer temperatures in marine ectotherms [71–74].

Regulation of gene expression

Temperature also had an evident effect on the regulation of gene expression in M. franciscanus gastrula embryos. GO terms identified as regulation of transcription by RNA polymerase II, translation regulator, and regulation of gene expression are relatively down-regulated in gastrula embryos raised under 17 °C relative to those raised under 13 °C (i.e., genes are comparatively up-regulated in the colder temperature treatment) (Fig. 3a). The turquoise WGCNA module, which was positively correlated to the 17 °C gastrula treatment and negatively correlated to the 13 °C gastrula treatment (Fig. 5), contained genes associated with transcription and several epigenetic processes. GO terms identified within the turquoise module included DNA templated transcription process, transcription coactivator, RNA polymerase binding, histone binding, chromatin binding, and RNA methylation (Fig. 6).

Epigenetic modifications, primarily consisting of DNA methylation, chromatin organization (e.g., histone and posttranslational modifications), and noncoding RNAs, are mechanisms of nongenetic variation by which the phenotypes of organisms can be altered faster than, and without the need of, changes in genotype [75, 76]. The dynamic changes in different epigenetic marks can act to regulate gene function without altering the DNA sequence, promoting phenotypic plasticity and potentially modulating the response to different environmental conditions [77–79]. Histone variants and modifications may activate or repress transcription processes by altering chromatin structures, impacting the regions of the genome that are available for transcription [80], and have been shown to mediate responses to changing environmental conditions in marine organisms [81–83]. Here, histone binding and chromatin remodeling genes are positively correlated with the 17 °C gastrula treatment. However, additional analyses such as ChIP-seq (i.e., chromatin immunoprecipitation sequencing) to locate regions targeted by histone modifications and DNA-binding proteins [84, 85] or ATAC-seq (i.e., assay for transposase-accessible chromatin with high-throughput sequencing) to assess genome-wide chromatin accessibility [86] are required to profile specific histone variants or modifications and their impact on gene expression.

Another layer to gene expression regulation includes the posttranscriptional modification of RNA molecules, an emerging concept known as RNA epigenetics or epitranscriptomics [87, 88]. The turquoise WGCNA module that was positively correlated to the 17 °C gastrula treatment included genes related to RNA modification, RNA methyltransferase, RNA methylation, and pseudouridine synthesis (Fig. 6). RNA modifications can occur in ribosomal RNA (rRNA), transfer RNA (tRNA), messenger RNA (mRNA) and small nuclear RNA (snRNA) [89, 90]. In addition to regulating gene expression, RNA modifications can modulate RNA transport and degradation, alternative splicing, protein binding, and development, among other essential biological processes [90–93]. RNA modifications such as the methylation of cytosine or adenosine are dynamic and reversible, allowing for alteration in protein-RNA interactions and rapid responses to changes in environmental conditions. For example, N⁶-methyladenosine (m6A), a very common RNA modification in eukaryotes, appears to play a role in the mammalian temperature stress response by promoting translation initiation of heat shock response genes [94]. Another common modification is pseudouridine (ψ), the first modified nucleoside discovered in RNA [95]. An up-regulation of genes related to pseudouridine synthesis could indicate enhanced transcript stability under elevated temperatures, as was observed in yeast undergoing heat shock [96]. Despite over 100 types of RNA modifications having been recorded [97], epitranscriptomic studies in invertebrates and marine species remain limited. Nonetheless, recent technological advances such as m6A-seq [98, 99] and ψ-seq [96] have provided exciting future avenues to explore epitranscriptomics in non-model systems.

Although evidence of transcription regulation was observed, it is difficult to conclude if this led to an increase or decrease of gene expression, particularly with respect to the functional significance of histone and RNA modifications. These mechanisms are much less studied in marine invertebrates than other modifications such as DNA methylation [77, 78], and although our data support that these modifications occurred in response to the environment, additional approaches are required to determine the precise modifications, their locations, and their impact on gene expression. In this study, there were more up-regulated than down-regulated genes in the 17 °C versus the 13 °C gastrula treatment. Nevertheless, future studies pairing comparative epigenetic and epitranscriptomic analyses with transcriptomic approaches are required to elucidate how these various mechanisms influence gene expression in response to different environmental conditions in M. franciscanus.

p CO ₂ influenced gene expression and body size of gastrula embryos

The pCO₂ treatment influenced gene expression patterns at the gastrula stage, although to a lesser degree than temperature, and the interaction between the two factors was not significant (Fig. 1d). Gastrula pCO₂ conditions explained 12.5% of the observed variance (p = 0.034) with 49 genes up-regulated and 202 genes down-regulated in embryos raised under 1050 µatm pCO₂ relative to those raised under 475 µatm pCO₂ (Figs. 1 and 2). In this study, we anticipated that the 475 µatm pCO₂ treatment was not stressful, as it represented the average ambient pCO₂ levels M. franciscanus regularly experience in their natural habitat [29]. Evidence suggests that calcifying marine organisms such as M. franciscanus are sensitive to declines in ocean pH (i.e., increases in pCO₂ levels) [100–102], and while M. franciscanus may periodically experience elevated pCO₂ conditions in nature during upwelling events [15, 28, 30, 103], the 1050 µatm pCO₂ treatment was expected to induce a stress response.

While the effect of pCO₂ on gastrula gene expression patterns was less than the effect of temperature, there was a pronounced impact of pCO₂ conditions on gastrula body size. Gastrula raised under elevated pCO₂ conditions (i.e., 1050 µatm) were significantly smaller than those raised under 475 µatm [31]. Interestingly, WGCNA module-trait relationships showed that the trait for gastrula body size shared highly similar module correlations to that of the lower pCO₂ treatment (475 µatm) (Fig. 5). Therefore, pCO₂ elicited a relatively muted transcriptomic response compared to temperature, but appeared to have a much greater influence at the organismal level. Below we discuss the expression patterns of genes affected by the gastrula pCO₂ treatment, which included those related to metabolism, ion transport, and the cellular stress response.

Metabolism and ion transport

Gastrula raised under the lower pCO₂ treatment (i.e., 475 µatm) expressed genes associated with macromolecule biosynthetic processes. In contrast, those raised under the elevated pCO₂ treatment (i.e., 1050 µatm) expressed genes associated with macromolecule catabolic processes (Fig. 3b). This may, in part, explain the difference in body size that was observed as a result of the pCO₂ conditions. Gastrula in the 475 µatm pCO₂ treatment appeared to construct proteins and other macromolecules to maintain their growth and body size, while gastrula in the 1050 µatm pCO₂ treatment underwent catabolic processes, possibly to obtain the energy required to respond to pCO₂ stress. The pCO₂ stress response may include an increase of ion transport as a means of maintaining acid-base equilibrium given elevated H⁺ concentrations under high pCO₂ conditions. DE analyses and WGCNA revealed that gastrula in the 1050 µatm treatment up-regulated genes related to ion binding, active transmembrane transporter, and ATPase coupled to movement of substances (Figs. 3b and 6). Similarly, increased expression of ion transport genes has been observed in gastrula embryos of S. purpuratus exposed to moderately elevated pCO₂ levels (e.g., ~ 800 µatm) [104], although the transcriptomic response of S. purpuratus embryos to pCO₂ stress can be influenced by maternal effects [105].

Cellular stress response

Elevated pCO₂ exposure also often impacts expression of cellular stress response genes in marine metazoans [106]. The cellular stress response is associated with the increased synthesis of molecular chaperones [107], including heat shock proteins (HSPs). Here, gastrula embryos raised under 1050 µatm pCO₂ differentially up-regulated genes related to heat shock protein binding (Fig. 3b). This result conflicts with another study in M. franciscanus that found that expression of the molecular chaperone Hsp70 decreased under elevated pCO₂ levels [27], although this was measured at the larval stage. However, increased expression of HSP genes under elevated pCO₂ conditions has been observed in the Antarctic pteropod Limacina helicina antarctica [108], in juveniles of the stony coral Acropora millepora [109], and in the cold-water coral Desmophyllum dianthus [110]. HSPs act as molecular chaperones in the cellular stress response by assisting in protein transport, protein folding and unfolding, stabilization of denatured proteins, and degradation of misfolded proteins [111, 112]. Thus, gastrula in the 1050 µatm pCO₂ treatment appeared to respond to stressful conditions via increased expression of HSP genes to establish or maintain physiological equilibrium.

Temperature was the dominant factor at the prism stage

At the prism stage, temperature accounted for 22.8% of the observed variance (p < 0.001) with 4132 genes up-regulated and 4286 genes down-regulated in embryos raised under 17 °C relative to those raised under 13 °C (Figs. 1 and 2). Unlike at the gastrula stage, responses to temperature measured at the molecular level were also observable at the organismal level. Development at the warmer 17 °C treatment led to an increase in prism body size as well as a modest increase in prism thermotolerance [31]. Indeed, prism body size and thermotolerance traits exhibited highly similar WGCNA module correlations as the prism 17 °C treatment (Fig. 5). Therefore, the transcriptomic response to temperature appears to have influenced both growth and resistance to heat stress in M. franciscanus prism embryos.

Similarities to the transcriptomic response at the gastrula stage

Like the gastrula stage, prism embryos raised at 17 °C exhibited increased expression of genes related to oxidative stress, transmembrane transport, and metabolic processes. As previously discussed, oxidative stress caused by the production of ROS can lead to cellular damage of lipids, proteins, and nucleic acids [49, 50, 60]. GO enrichment from the DE analysis and WGCNA identified terms including oxidoreductase, antioxidant, response to oxidative stress, and response to reactive oxygen species from genes up-regulated in response to the warmer temperature treatment (Figs. 4 and 6). At the gastrula stage, embryos in the 17 °C treatment also expressed genes associated with protein and DNA repair, but there was no evidence of increased macromolecule repair gene expression at the prism stage. This could indicate that while both stages responded to oxidative stress due to warmer temperatures, there was less cellular damage of proteins and nucleic acids incurred by the prism stage.

Also similar to the gastrula stage, warmer temperature conditions caused an increased expression of genes related to transmembrane transport in prism embryos. Specifically, up-regulated genes were related to ion binding, active transmembrane transporter, ATP hydrolysis coupled cation transmembrane transport, and sodium and potassium ion transmembrane transporters (Figs. 4 and 6), which may indicate osmoregulation and the maintenance of homeostasis. Prism embryos in the 17 °C treatment also increased expression of genes related to energetic processes (e.g., ATP metabolic process, organic acid metabolic process, lipid catabolic process, cyclic nucleotide metabolic process, and carbohydrate metabolic process) possibly to generate the energy required to support active transmembrane transport of ions and other materials. The up-regulation of genes related to energy production may have also supported the increased growth of prism embryos under warmer temperatures. Genes within the blue WGCNA module, which was positively correlated to the 17 °C prism treatment, were related to amide and nucleoside phosphate biosynthetic processes. This may indicate anabolic metabolism in which additional synthesis of proteins and nucleic acids supported larger body sizes. In contrast, DE analysis revealed an increase in expression of genes related to the negative regulation of biological and metabolic processes in prism embryos raised at 13 °C. This supports the predicted expectation that organisms exhibit decreased metabolism under colder temperatures [70–74].

DNA damage and gene expression regulation at lower temperatures

In general, global change biology research in marine systems has focused on the negative consequences of increasing temperatures associated with ocean warming [113]. However, given the expected rise in variable and extreme weather events, the impact of decreased temperatures is also an important consideration, especially in regions dominated by upwelling. Both DE analysis and the WGCNA indicated that prism embryos in the 13 °C treatment increased expression of genes related to DNA damage and repair (Figs. 4 and 6). Specifically, identified GO terms included cellular response to DNA damage stimulus, DNA helicase, DNA metabolic process, DNA binding, and DNA polymerase. Increased DNA damage as a result of low temperature stress has been recorded in the Pacific white shrimp Litopenaeus vannamei [114]. Although 13 °C is within the range of temperatures that M. franciscanus experience in the Santa Barbara Channel where urchins were collected for this study, it is lower than the annual average of ~ 15 °C [29] and may have generated stress and cellular damage in the prism embryos, leading to the activation of repair mechanisms.

There also appeared to be a relative up-regulation of genes associated with gene expression regulation in prism embryos raised under the lower temperature treatment, including genes related to translation regulator, transcription coregulator, regulation of gene expression, and regulation of transcription by RNA polymerase II (Fig. 4). Additionally, the turquoise WGCNA module, which contained genes related to epigenetic and epitranscriptomic modifications, was positively correlated to the 13 °C prism treatment (Fig. 5). Gene expression in prism embryos raised at 13 °C appeared to be regulated by histone and chromatin modifications, which can influence genome accessibility, and by RNA modifications such as RNA base methylation via methyltransferase activity (Figs. 4 and 6). Interestingly, the turquoise WGCNA module was positively correlated to the warmer temperature treatment at the gastrula stage, but was positively correlated to the colder temperature treatment at the prism stage. Thus, how the M. franciscanus epigenome and epitrancriptome respond to temperature conditions can vary during development. This difference in gene expression regulation may have contributed to temperature affecting body size at one developmental stage and not at the other.

Temperature did not influence HSP gene expression

Higher levels of HSPs, which act as molecular chaperones that protect cells from stress-induced damage, have been shown to confer increased thermotolerance across a variety of marine taxa [115–118]. Therefore, we may have expected an up-regulation of HSP genes linked with the slight increase in thermotolerance measured at the prism stage [31]. However, genes related to HSPs were not differentially expressed as a result of temperature at either the gastrula or prism stages. Our results contrast with a study in S. purpuratus that found expression of Hsp70 and Hsp90 increased at higher temperatures [41]. It is important to note that our study implemented a comparative approach that did not detect a significant difference in the relative expression of HSP genes. This approach may not be capable of detecting more subtle changes in HSP expression, which can vary between genes. For example, in the Antarctic urchin Sterechinus neumayeri, thermal stress led to an increase in Hsp90 expression, while Hsp70 expression remained unchanged [119]. A quantitative approach (e.g., using qPCR) may provide more insight into the absolute expression of genes like Hsp70 and Hsp90 during M. franciscanus development.

In other investigations of sea urchin early development, increased Hsp70 expression was generally not observed under moderate warming scenarios. One study found that Hsp70 was not transcriptionally up-regulated in M. franciscanus until larvae were exposed to temperatures at or above 20 °C [27]. A study in S. purpuratus only found induction of Hsp70 occurred at temperatures above 21 °C [120]. Therefore, the 17 °C treatment may not have been extreme enough to induce differential expression of HSP genes. Furthermore, in the green sea urchin Psammechinus miliaris, expression of HSP genes were low during early development relative to expression in adults [121]. The authors suggested that HSP expression was limited during this time [121] because over-expression of HSPs could have negative consequences for successful early development [122]. Therefore, large increases in HSP expression may be restricted during M. franciscanus early development.

The prism stage exhibited a limited transcriptomic response to pCO₂

In other studies, echinoderms raised under elevated pCO₂ conditions have exhibited altered expression of genes related to skeletogenic pathways, spicule matrix proteins, cellular stress response, ion regulation and transport, apoptosis, metabolism and ATP production [104, 105, 123–126]. Here, the pCO₂ treatment had a relatively minimal effect on gene expression patterns, particularly at the prism stage. The pCO₂ treatment explained only 9.6% of the observed variance at this stage and was not significant (p = 0.084) (Fig. 1f). This contrasts with observations made at the organism level in which elevated pCO₂ resulted in smaller prism embryos, although this could be offset by the positive effect of temperature, which acted as the dominant factor influencing body size [31]. It is interesting that the transcriptomic response to elevated pCO₂ was more evident in gastrula than in prism embryos particularly because at the prism stage, skeletal rod formation occurs. We may expect increased stress associated with prism embryos undergoing calcification processes under lowered pH conditions, however, evidence of this was not detected by changes in gene expression.

It is possible that while there is a clear phenotypic difference in prism embryos raised under high versus low pCO₂ conditions, the transcriptomic changes underlying this difference are too subtle to be identified statistically. Alternatively, the prism stage may simply lack a robust transcriptional response to the 1050 µatm pCO₂ treatment. For instance, the Mediterranean sea urchin Paracentrotus lividus exhibits different transcriptomic responses depending on the magnitude of the pH stressor [127]. Decreased pH conditions caused P. lividus embryos to increase their expression of calcification genes, but not once the pH dropped below a certain threshold [127]. A similar result was observed in S. purpuratus in which embryos raised under a high pCO₂ treatment designed to reflect near-future levels was relatively muted relative to those raised under a more moderate pCO₂ treatment designed to reflect present-day low pH conditions [104]. The authors speculated that the transcriptional response required for acclimating to a more extreme pCO₂ level was too metabolically expensive, and the embryos instead opted to conserve energy to ensure short-term survival, perhaps until environmental conditions became more favorable [104]. While a failure of embryos to respond at the transcriptomic level may allow for continued successful development under high pCO₂ conditions, there may be important physiological consequences such as the observed reduction in body size [31]. Thus, the lack of a transcriptomic response to high pCO₂ may have important fitness consequences for M. franciscanus.

The genetic structure of the sea urchins used in this experiment may have contributed to the lack of an expected response to the pCO₂ treatment. In the green sea urchin Strongylocentrotus droebachiensis, a quantitative genetic breeding design implemented by Runcie and colleagues demonstrated that changes in gene expression as a result of differences in pH exposure were minor relative to gene expression differences as a result of parentage [40]. Thus, minimal transcriptomic responses to pCO₂ may be due to genetic variation in M. franciscanus. Furthermore, the environmental exposure history of the adult urchins may have generated non-genetic parental effects (i.e., transgenerational plasticity), which can also generate a limited transcriptomic response to high pCO₂ [105]. While all embryo cultures for this experiment were composed of the same mixture of progeny from a cross between one male and five females, it is possible that the sea urchins collected for this experiment may be from a population with a relatively muted transcriptomic response to high pCO₂ conditions.

It has been proposed that selection and local adaptation acts on populations that are regularly exposed to high pCO₂ conditions, such as those that often experience upwelling conditions within the California Current System (CCS) [15], and that these populations may harbor genotypes that are resistant to low pH conditions [124, 128–132]. In S. purpuratus, transcriptomic responses to high pCO₂ levels can vary by the frequency in which the sea urchin populations are exposed to upwelling conditions [132]. In particular, urchins from populations frequently exposed to low pH have greater transcriptomic responses to high pCO₂ than those that experience low pH less often [132]. The adult sea urchins used in this study were collected from a site in the Santa Barbara Channel. While this area does experience periods of low pH due to upwelling [28, 30], low pH events occur less frequently than at more northern sites within the CCS [15, 133, 134]. Therefore, the urchins used in this study may be comparatively less adapted towards mounting a transcriptomic response to high pCO₂.

Performance under current and future ocean conditions

Moderate ocean warming may be favorable for M. franciscanus early development by providing larger body sizes and increased thermotolerance at the prism stage [31]. The warmer temperature treatment could even mitigate the stunting effect of elevated pCO₂ on prism body size [31]. This effect of warmer temperatures, however, may only be beneficial on a short-term basis. Gene expression analyses indicated that embryos raised under 17 °C responded to cellular stress, and while there were no indications of negative impacts at the phenotypic level, there may be trade-offs and consequences to developing under warmer temperatures such as increased incidences of disease [135]. Prolonged heat exposure may eventually become detrimental, and negative carry-over effects can arise at later life stages [46, 47]. Additionally, the observed plasticity at 17 °C may not extend to more severe warming scenarios. For example, a study in adult M. franciscanus found that although mortality did not vary between urchins acclimated to 15 °C or 18 °C, mortality was significantly higher at a more extreme temperature of 21 °C [26]. Nevertheless, our study revealed that M. franciscanus may be quite resilient to warming, and may even benefit from relatively brief and modest increases in temperature during early development. At the urchin collection site, temperatures of 17 °C are currently recorded during the summer months [29], and in the future, this temperature is likely to be reached more often given unmitigated climate change. More research is required to determine how M. franciscanus will be impacted as ocean warming continues, particularly as marine heatwaves increase in frequency [136, 137].

During early development, M. franciscanus appear to be more susceptible to rising pCO₂ levels than to rising temperatures. The lack of a robust transcriptional response paired with a decrease in body size indicates that exposure to elevated pCO₂ is detrimental to developing embryos. Continued ocean acidification is therefore likely to have adverse impacts on future M. franciscanus populations, although this could be offset somewhat by the positive effects of simultaneous ocean warming [31, 138–140]. However, during seasonal upwelling events, M. franciscanus are exposed to corrosive pH conditions that lack the mitigating effects of warmer temperatures [15, 141]. The upwelling season in this region typically extends from early spring until late summer or fall, and is characterized by variable fluctuations between periods of upwelling and the relaxation of upwelling [141, 142]. This overlaps with the natural spawning period of M. franciscanus that occurs annually during spring and early summer months [143–145]. As such, spawning that occurs during or immediately prior to an upwelling event will subject developing embryos to high pCO₂ conditions paired with colder temperatures. Given ocean acidification and the increase in upwelling frequency and intensity that is predicted with continued climate change [146–148], the likelihood of M. franciscanus developing under stressful pCO₂ conditions should rise in the future.

The early developmental stages of M. franciscanus exhibited plasticity at both a transcriptomic and phenotypic level under different temperature and pCO₂ conditions. The extent of this plasticity, however, varied by developmental stage and by stressor. Although higher temperatures appeared to induce cellular stress, M. franciscanus exhibited a robust transcriptional response to temperature, and may even benefit from ocean warming. The transcriptomic response to pCO₂ is much more limited, and may therefore increase the vulnerability of this species to ocean acidification. Developing embryos of M. franciscanus are likely already negatively impacted by present-day low pH conditions that occur as a result of seasonal upwelling.

Molecular-level approaches such as the one used here are valuable means by which to assess how an organism is responding to its environment, particularly when paired with other trait data. For instance, if examining differences in body size alone, one may conclude that there was no effect of temperature on the gastrula stage and that they are either resilient or are unable to respond to the different temperatures used in this experiment. However, insight into gastrula transcriptomic patterns revealed that the embryos responded to temperature on a molecular level, possibly allowing them to maintain their growth and resulting body size under either temperature treatment. This underscores the advantages of using an integrative approach that spans multiple hierarchical levels (i.e., from molecules to organisms to ecosystems) whenever possible to obtain a more comprehensive representation of how changing environments affect living systems.

Accurately predicting how organisms will respond to future ocean conditions is necessary for the implementation of effective conservation and fisheries management strategies. This study provides essential insight towards understanding how an ecologically valuable fishery species is affected by climate change. Although we found M. franciscanus exhibits plasticity under temperature and pCO₂ stress, more work is required to determine how environmental stressors will impact the overall fitness and adaptive potential of natural populations. Detrimental impacts that occur during early development can reduce the quantity or quality of juvenile recruits. Poor recruitment will likely not be immediately evident to the M. franciscanus fishery until there is a marked decrease in the number of mature, harvestable adults. This delay in noticeable population declines may cause substantial alterations to coastal macroalgae ecosystems and financial losses to the fishing industry. Studies such as the one presented here are vital for developing proactive and adaptive management strategies to ensure climate-ready fisheries [149, 150].

Animal collection and culturing

Red sea urchins were collected and spawned as described in [31]. Briefly, adults were collected from Ellwood Mesa, Goleta, California, USA (34° 25.065’N, 119° 54.092’W) at 14-m depth via SCUBA on February 21, 2018 under California Scientific Collecting Permit SC-1223 and transported to the Marine Science Institute at the University of California Santa Barbara (UCSB). Spawning was induced by injecting 0.53 M KCl into the coelom through the perioral membrane [145]. Eggs from five individual females and sperm from a single male were collected. A subsample of eggs from each female was fertilized with sperm from the male and high fertilization success was examined for each cross (i.e., visually confirming the formation of fertilization envelopes). These subsamples were only used to verify suitable male-female compatibility and were discarded prior to the experiment. An approximately equal number of eggs from each of the five females were gently pooled together. The pool of eggs was fertilized by slowly adding dilute, activated sperm from the male until approximately 98% fertilization success was reached. Performing crosses with a single male ensured that all cultures were composed of full- or half-sibling embryos. This approach was selected in an effort to limit paternal genetic variability and differences in male-female interactions that could otherwise impact the results of the study.

As described in [31], the newly-fertilized embryos were raised in three replicate culture vessels for each of four different combined temperature and pCO₂ treatments: 1) 17 °C and 1050 µatm pCO₂, 2) 17 °C and 475 µatm pCO₂, 3) 13 °C and 1050 µatm pCO₂, and 4) 13 °C and 475 µatm pCO₂. These treatments represent current conditions measured in the Santa Barbara Channel (SBC) as well as projected future ocean conditions given continuing warming and acidification. It should be noted that the SBC is within a highly dynamic region in which water conditions can vary greatly. For example, the temperature recorded where and when the M. franiscanus were collected for this study was 13.3 °C, but throughout the year prior, the temperature at the collection site fluctuated between 10.0–20.5 °C, with an average temperature of 15.1 °C [29]. Static treatments were selected, however, because the early development of M. franciscanus occurs rapidly (e.g., < 2.5 days in this study) so unlike their adult counterparts, the embryos are likely to experience relatively stable conditions. The treatment of 13 °C and 475 µatm pCO₂ represents conditions that are regularly measured in the SBC [29]. 13 °C is on the lower end of the range of temperatures observed in this region, but is regularly observed during seasonal upwelling events. During upwelling events, combinations of decreased temperatures and elevated pCO₂ levels are common [15, 28, 30, 103]. These upwelling conditions are represented by the combined 13 °C and 1050 µatm pCO₂ treatment. Water conditions in this area also currently reach elevated temperatures of 17 °C on occasion [29], but as ocean warming continues, these conditions are expected to increase in frequency [151]. The combined 17 °C and 1050 µatm pCO₂ treatment represents future conditions as both ocean warming and ocean acidification progress. The temperature, pH, salinity, and total alkalinity were measured in each culture vessel daily as detailed in [31].

Sample generation and sequencing

The early gastrula stage was designated by the formation of secondary mesenchyme cells and the extension of the archenteron to approximately one-half the embryo’s body length (~ 23.5 hours post-fertilization (hpf) at 17 °C and ~ 32.5 hpf at 13 °C). The prism stage was designated by the tripartitioning of the archenteron, the formation of skeletal rods, and the pyramid-like body shape of the embryo (~ 44 hpf at 17 °C and ~ 55.5 hpf at 13 °C). Embryos at both developmental stages were sampled from each culture vessel by gently concentrating the embryos onto a submerged 35-µm mesh filter and transferring them into a 15-mL Falcon tube using a plastic transfer pipet. The concentration of embryos per mL of FSW was calculated by counting three aliquots of embryos so that a coefficient of variance (CV) of less than 10% was reached. At both the gastrula and prism stages, 5000 embryos from each culture vessel were transferred into a 1.5 mL microcentrifuge tube. Embryos were quickly pelleted by centrifugation, excess seawater was removed, and the samples were flash frozen in liquid nitrogen. The samples were stored at -80 °C until RNA extractions were performed.

Total RNA extractions were performed on samples from each culture vessel at both the gastrula and prism stages for a total of 24 RNA extractions. RNA extractions were performed by adding 500 µL of Trizol^® reagent (Invitrogen) to each sample, and passing the entire contents through needles of decreasing sizes (i.e., three passes through a 21-gauge needle, a 23-gauge needle, and then a 25-gauge needle) until homogenized. The RNA was isolated using a chloroform addition. The RNA was precipitated in 100% isopropyl alcohol, washed in ice cold 80% ethanol, and resuspended in DEPC-treated water. RNA purity, quantity, and quality were verified using a NanoDrop^→ ND100, a Qubit^Ⓡ fluorometer, and a TapeStation 2200 system (Agilent Technologies). The RNA samples were submitted to the Genome Center at the University of California, Davis where 24 libraries were generated using poly(A) enrichment. The 24 libraries were pooled and sequenced across two lanes on an Illumina HiSeq 4000 with 100 base-pair (bp) single reads.

Data processing

Adapter sequence contamination and base pairs with quality scores below 30 were removed from the raw sequence data using Trim Galore! (version 0.4.1) [152]. FastQC (version 0.11.5) [32] was used to verify sequence quality. The trimmed sequence data were mapped onto a developmental transcriptome for M. franciscanus [33] and expression values were calculated using RSEM (version 1.3.0) [34] and bowtie2 (version 2.3.2) [153]. The limma package [35] in R (version 3.4.4) was used to filter the sequence data to those that had at least 0.5 counts per million mapped reads across at least three of the samples. A scale normalization was applied to the read counts using a trimmed mean of M-values (TMM) normalization method [154]. limma was used to voom-transform the data, converting the read counts to log-counts per million while accounting for sample-specific quality weights and the blocking design of the experiment. Principal component analyses (PCAs) were performed on the filtered, normalized, and voom-transformed data to examine distances between the samples using the prcomp function in R. Permutational multivariate ANOVAs were performed using the adonis function in the package vegan [155] to establish the proportion of variance explained by fixed factors (i.e., developmental stage, temperature treatment, and pCO₂ treatment).

DE Analysis and WGCNA

Because gene expression patterns in M. franciscanus vary significantly by developmental stage [33] and because comparing gene expression patterns across development was not a goal of the current study, gastrula stage and prism stage DE analyses were executed separately. DE analyses were performed in limma by making pairwise comparisons between temperature treatments (17 °C versus 13 °C) and between pCO₂ treatments (1050 µatm versus 475 µatm) for each developmental stage. Moderated t-statistics and p-values adjusted by the Benjamini and Hochberg’s method were used to control the false discovery rate [156]. Differentially expressed genes (genes relatively down- or up-regulated) were determined for each pairwise comparison of interest (adjusted p-value < 0.05). Gene ontology (GO) analyses were performed following the GO_MWU package [36] in R. This package uses a Mann-Whitney U test to determine whether GO categories are significantly enriched by differentially regulated genes (up- or down-regulated). For each DE analysis, values of log-fold change in expression for all genes were used to test GO enrichment within three categories: molecular function (MF), biological process (BP), and cellular component (CC).

WGCNA [37] was performed in R to identify clusters of similarly expressed genes which were sorted into modules (minimum of 30 genes per module). A signed network was constructed and used to generate modules with highly correlated eigengenes that were merged using a threshold of 0.26. Eigengene expression was then correlated with each trait. Traits were defined as gastrula temperature (17 or 13 °C), gastrula pCO₂ (1050 or 475 µatm), prism temperature (17 or 13 °C), and prism pCO₂ (1050 or 475 µatm). Because the experimental treatments have been shown to affect embryo body size and thermotolerance [31], these phenotypic traits were also incorporated in the WGCNA analysis. Gastrula and prism body size traits were included using average length values for each of the 24 samples. Thermotolerance, which was only measured for the prism stage, was included using the average LT₅₀ value of each treatment, which corresponds to the temperature at which 50% mortality occurs. The methods of how body size and thermotolerance values were obtained are described in detail in [31]. A heatmap was generated to visualize significant correlations (-1 to 1, p-value < 0.05) between each trait and module eigengene. GO analysis on WGCNA modules was performed following the GO_MWU package [36] using module membership values (i.e., kME values) for all genes across gastrula and prism stages. GO terms within MF, BP, and CC categories were identified for each module eigengene.

pCO₂

Partial pressure of carbon dioxide; RNA-seq:RNA sequencing; PCA:Principal component analysis; PC1:Principal component 1; PC2:Principal component 2; DE:Differential expression; GO:Gene ontology; MF:Molecular function; BP:Biological process; CC:Cellular component; WGCNA:Weighted gene co-expression network analysis; mm:Millimeters; ROS:Reactive oxygen species; ChIP-seq:Chromatin immunoprecipitation sequencing; ATAC-seq:Assay for transposase-accessible chromatin with high-throughput sequencing; rRNA:Ribosomal RNA; tRNA:Transfer RNA; mRNA:Messenger RNA; snRNA:Small nuclear RNA; m6A:N⁶-methyladenosine;ψ:Pseudouridine; HSP:Heat shock protein; CCS:California Current System; UCSB:University of California Santa Barbara; SBC:Santa Barbara Channel; hpf:Hours post-fertilization; CV:Coefficient of variance; bp:Base-pair; TMM:Trimmed mean of M-values

Ethics approval and consent to participate

This study was conducted on a marine invertebrate that does not require Institutional Animal Care and Use Committee (IACU) protocols or approval. However, the specimens were collected and handled in accordance with the requirements of a California Scientific Collecting Permit granted to GEH (SC-1223) that has been noted in the Methods. In addition, no human research subjects were used for this study, and IRB approval was not required.

Consent for publication

Not applicable.

Availability of data and material

The datasets generated and/or analyzed during the current study are available in the NCBI Short Read Archive under Bioproject accession number PRJNA637102. Additional R data and analysis scripts are available at a Github repository (https://github.com/julietmwong27/Mfranciscanus_RNAseq_Ranalysis).

Competing interests

The authors declare that they have no competing interests.

Funding

This research was supported by funds from the UC Climate Champion award from the University of California to GEH. The funding source did not influence the design or execution of the study and the authors acted independently throughout the course of the research project.

Animal collections were supported by resources from the Santa Barbara Coastal Long-Term Ecological Research program (NSF award OCE-1232779). During this project, JMW was supported by a U.S. National Science Foundation (NSF) Graduate Research Fellowship under Grant No. 1650114. Analytical and writing phases of the project were supported by NSF award IOS-1656262 to GEH. The authors acknowledge the use of the UCSB and UCOP-supported Biological Nanostructures Laboratory within the California NanoSystems Institute. This work also used computing resources supported by NSF under grant No. ABI-1458641 to Indiana University.

Authors’ contributions

JMW designed and implemented the experiment. Collection of the samples, RNA isolation, and bioinformatic analyses were performed by JMW. JMW and GEH analyzed the results and wrote the manuscript. Both authors have read and approved the final version.

Acknowledgements

We thank Christoph Pierre, Director of Marine Operations at UC Santa Barbara, for his assistance with boating and sea urchin collections. We would also like to thank Logan Kozal, Terence Leach, Maddie Housh, Jannine Chamorro, Dr. Marie Strader, and Dr. Umihiko Hoshijima for their assistance during the experiment.

Ebert TA, Dixon JD, Schroeter SC, Kalvass PE, Richmond NT, Woodby WA. DA. Growth and mortality of red sea urchins Strongylocentrotus franciscanus across a latitudinal gradient. Mar Ecol Prog Ser. 1999;190:189–209.
Rogers-Bennett L. The ecology of Strongylocentrotus franciscanus and Strongylocentrotus purpuratus. Developments in Aquaculture Fisheries Science vol. 2007;37:393–425.
Leighton D, Jones L, North W. Ecological relationships between giant kelp and sea urchins in southern California. In: Proceedings of the 5th international seaweed symposium. Edited by Young E, Maclachlan J. Pergamon Press. 1966. 141–153.
Filbee-Dexter K, Scheibling RE. Sea urchin barrens as alternative stable states of collapsed kelp ecosystems. Mar Ecol Prog Ser. 2014;495:1–25. doi:10.3354/meps10573.
Tegner MJ, Levin LA. Spiny lobsters and sea urchins: Analysis of a predator-prey interaction. J Exp Mar Biol Ecol. 1983;73(2):125–50.
McLean JH. Sublittoral ecology of kelp beds of the open coast area near Carmel, California. Biol Bull. 1962;122(1):95–114.
A EJ FPJ. Sea otters: their role in structuring nearshore communities. Science. 1974, 185(4156):1058–1060.
Andrew N, Agatsuma Y, Ballesteros E, Bazhin A, Creaser E, Barnes D, Botsford L, Bradbury A, Campbell AK, Dixon J, et al. Status and management of world sea urchin fisheries. Oceanogr Mar Biol: Annu Rev. 2002;40:343–425.
Keesing J, Hall K. Review of harvests and status of world sea urchin fisheries points to opportunities for aquaculture. J Shellfish Res. 1998;17(5):1597–604.
The Pacific Fisheries Information Network (PacFIN). https://pacfin.psmfc.org. Accessed 29 May 2020.
Rogers-Bennett L, Okamoto D. Mesocentrotus franciscanus and Strongylocentrotus purpuratus. In: Developments in Aquaculture and Fisheries Science. vol. 43. 2020. 593–608.
Sato KN, Powell J, Rudie D, Levin LA. Evaluating the promise and pitfalls of a potential climate change–tolerant sea urchin fishery in southern California. ICES J Mar Sci. 2018, 75(3):1029–1041. doi:10.1093/icesjms/fsx225.
Gentemann CL, Fewings MR, García-Reyes M. Satellite sea surface temperatures along the West Coast of the United States during the 2014–2016 northeast Pacific marine heat wave. Geophys Res Lett. 2017;44:312–9. doi:10.1002/2016GL071039.
Oliver ECJ, Donat MG, Burrows MT, Moore PJ, Smale DA, Alexander LV, Benthuysen JA, Feng M, Gupta AS, Hobday AJ, et al. Longer and more frequent marine heatwaves over the past century. Nat Commun. 2018, 9(1324). doi:10.1038/s41467-018-03732-9.
Chan F, Barth JA, Blanchette C, Byrne RH, Chavez F, Cheriton O, Feely RA, Friederich G, Gaylord B, Gouhier T, et al. Persistent spatial structuring of coastal ocean acidification in the California Current System. Scientific Reports. 2017, 7(2526). doi:https://doi.org/10.1038/s41598-017-02777-y.
Daigle RM, Metaxas A. Vertical distribution of marine invertebrate larvae in response to thermal stratification in the laboratory. J Exp Mar Biol Ecol. 2011;409:89–98.
Maboloc EA, Batzel G, Grünbaum D, Chan KYK. Vertical distribution of echinoid larvae in pH stratified water columns. Mar Biol. 2020, 167.
Byrne M. Impact of ocean warming and ocean acidification on marine invertebrate life history stages: Vulnerabilities and potential for persistence in a changing ocean. Oceanogr Mar Biol: Annu Rev. 2011;49:1–42.
Dupont S, Thorndyke M. Impact of CO₂-driven ocean acidification on invertebrates early life-history – What we know, what we need to know and what we can do. Biogeosciences. 2009;6:3109–31.
Gosselin L, Qian P-Y. Juvenile mortality in benthic marine invertebrates. Mar Ecol Prog Ser. 1997;146:265–82.
Kurihara H. Effects of CO₂-driven ocean acidification on the early developmental stages of invertebrates. Mar Ecol Prog Ser. 2008;373:275–84.
Byrne M. Global change ecotoxicology: Identification of early life history bottlenecks in marine invertebrates, variable species responses and variable experimental approaches. Mar Environ Res. 2012;76:3–15.
Byrne M, Przeslawski R. Multistressor impacts of warming and acidification of the ocean on marine invertebrates’ life histories. Integr Comp Biol. 2013;53(4):582–96. doi:10.1093/icb/ict049.
Frieder CA. Present-day nearshore pH differentially depresses fertilization in congeneric sea urchins. Biol Bull. 2014;226:1–7.
Reuter KE, Lotterhos KE, Crim RN, Thompson CA, Harley CDG. Elevated pCO₂ increases sperm limitation and risk of polyspermy in the red sea urchin Strongylocentrotus franciscanus. Glob Chang Biol. 2011, 17:163–171. doi:10.1111/j.1365-2486.2010.02216.x.
Hernández M, Bückle F, Guisado C, Barón B, Estavillo N. Critical thermal maximum and osmotic pressure of the red sea urchin Strongylocentrotus franciscanus acclimated at different temperatures. J Therm Biol. 2004;29:231–6.
O'Donnell M, Hammond L, Hofmann G. Predicted impact of ocean acidification on a marine invertebrate: elevated CO₂ alters response to thermal stress in sea urchin larvae. Mar Biol. 2009;156:439–46.
Kapsenberg L, Hofmann GE. Ocean pH time-series and drivers of variability along the northern Channel Islands, California, USA. Limnol Oceanogr. 2016;61:953–68. doi:10.1002/lno.10264.
SBC LTER: Reef: Bottom Temperature: Continuous water temperature, ongoing since 2000 ver 24. Environmental Data Initiative

https://doi.org/10.6073/pasta/5f8f6d12d9c8df2c26d6061bc2972599

Reed D. SBC LTER: Reef: Bottom Temperature: Continuous water temperature, ongoing since 2000 ver 24. Environmental Data Initiative. 2020. https://doi.org/10.6073/pasta/5f8f6d12d9c8df2c26d6061bc2972599.
Rivest EB, O'Brien M, Kapsenberg L, Gotschalk CC, Blanchette CA, Hoshijima U, Hofmann GE. Beyond the benchtop and the benthos: Dataset management planning and design for time series of ocean carbonate chemistry associated with Durafet®-based pH sensors. Ecol Inform. 2016;36:209–20. doi:https://doi.org/10.1016/j.ecoinf.2016.08.005.
Wong JM, Hofmann GE. The effects of temperature and pCO₂ on the size, thermal tolerance and metabolic rate of the red sea urchin (Mesocentrotus franciscanus) during early development. Mar Biol 2020, 167(33). doi:10.1007/s00227-019-3633-y.
Andrews S. FASTQC: a quality control tool for high throughput sequence data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc. 2010.
Wong JM, Gaitán-Espitia JD, Hofmann GE. Transcriptional profiles of early stage red sea urchins (Mesocentrotus franciscanus) reveal differential regulation of gene expression across development. Marine Genomics. 2019, 48(100692). doi:10.1016/j.margen.2019.05.007.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. Bmc Bioinformatics. 2011, 12(323).
Ritchie M, Phipson B, Wu D, Hu Y, Law C, Shi W, Smyth G. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Wright RM, Aglyamova G, Meyer E, Matz MV. Gene expression associated with white syndromes in a reef-building coral, Acropora hyacinthus. Bmc Genomics. 2015, 16(371). doi:10.1186/s12864-015-1540-2.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformatics 2008, 9(559). doi:10.1186/1471-2105-9-559.
Stumpp M, Wren J, Melzner F, Thorndyke M, Dupont S. CO₂ induced seawater acidification impacts sea urchin larval development I: Elevated metabolic rates decrease scope for growth and induce developmental delay. Comp Biochem Physiol A Mol Integr Physiol. 2011;160:331–40.
Fujisawa H. Differences in temperature dependence of early development of sea urchins with different growing seasons. Biol Bull. 1989;176(2):96–102.
Runcie D, Dorey N, Garfield D, Stumpp M, Dupont S, Wray G. Genomic characterization of the evolutionary potential of the sea urchin Strongylocentrotus droebachiensis facing ocean acidification. Genome Biol Evol. 2016;8(12):3672–84. doi:10.1093/gbe/evw272.
Runcie D, Garfield D, Babbitt C, Wygoda J, Mukherjee S, Wray G. Genetics of gene expression responses to temperature stress in a sea urchin gene network. Mol Ecol. 2012, 21(18). doi:10.1111/j.1365-294X.2012.05717.x.
Dan K, Okazaki K. Cyto-embryological studies of sea urchins. III. Role of the secondary mesenchyme cells in the formation of the primitive gut in sea urchin larvae. Biol Bull. 1956;110(1):29–42.
Wolpert L. Gastrulation and the evolution of development. Development. 1992;116(Suppl):7–13.
Burke RD. Morphogenesis of the digestive tract of the pluteus larva of Strongylocentrotus purpuratus: shaping and bending. Int J Invertebr Repr. 1980;2(1):13–21. doi:10.1080/01651269.1980.10553338.
Ettensohn CA, Malinda KM. Size regulation and morphogenesis: a cellular analysis of skeletogenesis in the sea urchin embryo. Development. 1993;119:155–67.
Beckerman A, Benton T, Ranta E, Kaitala V, Lundberg P. Population dynamic consequences of delayed life-history effects. Trends Ecol Evol. 2002;17(6):263–9.
Hettinger A, Sanford E, Hill TM, Russell AD, Sato KN, Hoey J, Forsch M, Page HN, Gaylord B. Persistent carry-over effects of planktonic exposure to ocean acidification in the Olympia oyster. Ecology. 2012;93:2758–68.
Seneca FO, Palumbi S. The role of transcriptome resilience in resistance of corals to bleaching. Mol Ecol. 2015;24:1467–84.
Ahmed RG. Is there a balance between oxidative stress and antioxidant defense system during develoopment? Med J Islamic World Acad Sci. 2005;15(2):55–63.
Costantini D. Oxidative stress in ecology and evolution: lessons from avian studies. Ecol Lett. 2008;11:1238–51. doi:10.1111/j.1461-0248.2008.01246.x.
Lesser MP. Elevated temperatures and ultraviolet radiation cause oxidative stress and inhibit photosynthesis in symbiotic dinoflagellates. Limnol Oceanogr. 1996;41(2):271–83.
Ali MB, Hahn E-J, Paek K-Y. Effects of temperature on oxidative stress defense systems, lipid peroxidation and lipoxygenase activity in Phalaenopsis. Plant Physiol Biochem. 2005;43:213–23. doi:10.1016/j.plaphy.2005.01.007.
Verlecar XN, Jena KB, Chainy GBN. Biochemical markers of oxidative stress in Perna viridis exposed to mercury and temperature. Chem Biol. 2007;167:219–26. doi:10.1016/j.cbi.2007.01.018.
Madeira D, Narciso L, Cabral HN, Vinagre C, Diniz MS. Influence of temperature in thermal and oxidative stress responses in estuarine fish. Comp Biochem Physiol A Mol Integr Physiol. 2013;116:237–43. doi:10.1016/j.cbpa.2013.06.008.
Lushchak VI, Bagnyukova TV. Temperature increase results in oxidative stress in goldfish tissues. 2. Antioxidant and associated enzymes. Comp Biochem Phys C. 2006;143:36–41. doi:10.1016/j.cbpc.2005.11.018.
Younus H. Oxidoreductases: Overview and Practical Application. In: Biocatalysis. Edited by Husain Q, Ullah M. Springer, Cham. 2019.
Genestra M. Oxyl radicals, redox-sensitive signalling cascades and antioxidants. Cell Signal. 2007;19(9):1807–19. doi:10.1016/j.cellsig.2007.04.009.
Esterházy D, King MS, Yakovlev G, Hirst J. Production of reactive oxygen species by complex I (NADH:ubiquinone oxidoreductase) from Escherichia coli and comparison to the enzyme from mitochondria. Biochemistry. 2008;47:3964–71.
Ngoka LCM. Dramatic down-regulation of oxidoreductases in human hepatocellular carcinoma hepG2 cells: proteomics and gene ontology unveiling new frontiers in cancer enzymology. Proteome Sci. 2008, 6(29). doi:10.1186/1477-5956-6-29.
Bandyopadhyay U, Das D, Banerjee RK. Reactive oxygen species: Oxidative damage and pathogenesis. Curr Sci. 1999;77(5):658–66.
Adams CJ, Kopp MC, Larburu N, Nowak PR, Ali MMU. Structure and molecular mechanism of ER stress signaling by the unfolded protein response signal activator IRE1. Front Mol Biosci. 2019, 6. doi:10.3389/fmolb.2019.00011.
Hetz C. The unfolded protein response: controlling cell fate decisions under ER stress and beyond. Nat Rev Mol Cell Biol. 2012;13:89–102.
Shental-Bechor D, Levy Y. Effect of glycosylation on protein folding: A close look at thermodynamic stabilization. PNAS. 2008;105:8256–61.
Lis H, Sharon N. Protein glycosylation: structural and functional aspects. Eur J Biochem. 1993;218:1–27.
DeSalvo MK, Sunagawa S, Voolstra CR, Medina M. Transcriptomic responses to heat stress and bleaching in the elkhorn coral Acropora palmata. Mar Ecol Prog Ser. 2010;402:97–113.
Li X, Heyer W-D. Homologous recombination in DNA repair and DNA damage tolerance. Cell Research. 2008;18:99–113.
Siudeja K, Bardin AJ. Somatic recombination in adult tissues: What is there to learn? Fly. 2017, 11(2):121–128. doi:10.1080/19336934.2016.1249073.
Vergara-Amado J, Silva AX, Manzi C, Nespolo RF, Cárdenas L. Differential expression of stress candidate genes for thermal tolerance in the sea urchin Loxechinus albus. J Therm Biol. 2017;68:104–9.
Leong PKK, Manahan DT. Metablic importance of Na⁺/K⁺-ATPase activity during sea urchin development. J Exp Biol. 1997;200:2881–92.
Gillooly JF, Brown JH, West GB, Savage VM, Charnov EL. Effects of size and temperature on metabolic rate. Science. 2001;293(5538):2248–51. doi:10.1126/science.1061967.
Arnberg M, Calosi P, Spicer J, Tandberg A, Nilsen M, Westerlund S, Bechmann R. Elevated temperature elicits greater effects than decreased pH on the development, feeding and metabolism of northern shrimp (Pandalus borealis) larvae. Mar Biol. 2013;160:2037–48. doi:10.1007/s00227-012-2072-9.
Manríquez PH, Torres R, Matson PG, Lee MR, Jara ME, Seguel ME, Sepúlveda F, Pereira L. Effects of ocean warming and acidification on the early benthic ontogeny of an ecologically and economically important echinoderm. Mar Ecol Prog Ser. 2017;563:169–84. doi:10.3354/meps11973.
McElroy DJ, Nguyen HD, Byrne M. Respiratory response of the intertidal seastar Parvulastra exigua to contemporary and near-future pulses of warming and hypercapnia. J Exp Mar Biol Ecol. 2012;416–417:1–7. doi:10.1016/j.jembe.2012.02.003.
Pimentel MS, Trübenbach K, Faleiro F, Boavida-Portugal J, Repolho T, Rosa R. Impact of ocean warming on the early ontogeny of cephalopods: a metabolic approach. Mar Biol. 2012;159:2051–9. doi:10.1007/s00227-012-1991-9.
Bondurianski R, Day T. Nongenetic inheritance and its evolutionary implications. Annu Rev Ecol Syst. 2009;40:103–25. doi:10.1146/annurev.ecolsys.39.110707.173441.
Deans C, Maggert KA. What do you mean. “Epigenetic”? Genetics. 2015;199:887–96. doi:10.1534/genetics.114.173492.
Eirin-Lopez JM, Putnam HM. Marine environmental epigenetics. Annu Rev Mar Sci. 2019;11:7.1–34. doi:10.1146/annurev-marine-010318-095114.
Hofmann GE. Ecological epigenetics in marine metazoans. Front Mar Sci. 2017, 4(4). doi:10.3389/fmars.2017.00004.
Verhoeven KJ, Vonholdt BM, Sork VL. Epigenetics in ecology and evolution: what we know and what we need to know. Mol Ecol. 2016;25:1631–8. doi:10.1111/mec.13617.
Berger SL. The complex language of chromatin regulation during transcription. Nature. 2007;447:407–12.
Gonzalez-Romero R, Suarez-Ulloa V, Rodriguez-Casariego J, Garcia-Souto D, Diaz G, Smith A, Pasantes JJ, Rand G, Eirin-Lopez JM. Effects of Florida Red Tides on histone variant expression and DNA methylation in the Eastern oyster Crassostrea virginica. Aquat Toxicol. 2017;186:196–204.
Rodriguez-Casariego JA, Ladd MC, Shantz AA, Lopes C, Cheema MS, Kim B, Roberts SB, Fourqurean JW, Ausio J, Burkepile DE, et al. Coral epigenetic responses to nutrient stress: Histone H2A.X phosphorylation dynamics and DNA methylation in the staghorn coral Acropora cervicornis. Ecol Evol. 2018;8(23):12193–207.
Veluchamy A, Rastogi A, Lin X, Lombard B, Murik O, Thomas Y, Dingli F, Rivarola M, Ott S, Liu X, et al. An integrative analysis of post-translational histone modifications in the marine diatom Phaeodactylum tricornutum. Genome Biol 2015, 16(102). doi:10.1186/s13059-015-0671-8.
Park PJ. ChIP–seq: advantages and challenges of a maturing technology. Nat Rev Genet. 2009;10:669–80.
O’Geen H, Echipare L, Farnham PJ. Using ChIP-Seq technology to generate high-resolution profiles of histone modifications. In: Epigenetics Protocols. Edited by Tollefsbol T, vol. 791. Humana Press. 2011. 265–286.
Buenrostro JD, Wu B, Chang HY, Greenleaf WJ. ATAC-seq: A method for assaying chromatin accessibility genome-wide. Curr Protoc Mol Biol 2015, 109:21.29.21–21.29.29. doi:10.1002/0471142727.mb2129s109.
He CRNA. epigenetics? Nat Chem Biol. 2000;6:863–5.
Saletore Y, Meyer K, Korlach J, Vilfan ID, Jaffrey S, Mason CE. The birth of the Epitranscriptome: deciphering the function of RNA modifications. Genome Biol. 2012, 13.
Romano G, Veneziano D, Nigita G, Nana-Sinkam SP. RNA methylation in ncRNA: Classes, detection, and molecular associations. Front Genet. 9. doi:10.3389/fgene.2018.00243.
Li S, Mason CE. The pivotal regulatory landscape of RNA modifications. Annu Rev Genom Hum G. 2014;15:127–50. doi:10.1146/annurev-genom-090413-025405.
Peer E, Rechavi G, Dominissini D. Epitranscriptomics: regulation of mRNA metabolism through modifications. Curr Opin Chem Biol. 2017;41:93–8. doi:10.1016/j.cbpa.2017.10.008.
Jia G, Fu Y, He C. Reversible RNA adenosine methylation in biological regulation. Trends Genet. 2013, 29.
Lewis CJT, Pan T, Kalsotra A. RNA modifications and structures cooperate to guide RNA–protein interactions. Nat Rev Mol Cell Biol. 2017;18:202–10.
Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian S-B. Dynamic m6A mRNA methylation directs translational control of heat shock response. Nature. 2015, 526:591–594. doi:10.1038/nature15377.
Charette M, Gray MW. Pseudouridine in RNA: What, where, how, and why. IUBMB Life. 2000;49:341–51.
Schwartz S, Bernstein DA, Mumbach MR, Jovanovic M, Herbst RH, León-Ricardo BX, Engreitz JM, Guttman M, Satija R, Lander ES, et al. Transcriptome-wide mapping reveals widespread dynamic-regulated pseudouridylation of ncRNA and mRNA. Cell. 2014;159:148–62.
Cantara WA, Crain PF, Rozenski J, McCloskey JA, Harris KA, Zhang X, Vendeix FAP, Fabris D, Agris PF. The RNA modification database, RNAMDB: 2011 update. Nucleic Acids Res. 2011;39:D195–201. doi:10.1093/nar/gkq1028.
Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, Cesarkas K, Jacob-Hirsch J, Amariglio N, Kupiec M, et al. Topology of the human and mouse m⁶A RNA methylomes revealed by m⁶A-sEq. Nature. 2012;485:201–6.
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons. Cell. 2012;149:1635–46.
Doney SC, Fabry VJ, Feely RA, Kleypas JA. Ocean acidification: The other CO₂ problem. Annu Rev Mar Sci. 2009;1:169–92. doi:10.1146/annurev.marine.010908.163834.
Hofmann GE, Barry JP, Edmunds PJ, Gates RD, Hutchins DA, Klinger T, Sewell MA. The effect of ocean acidification on calcifying organisms in marine ecosystems: An organism-to-ecosystem perspective. Annu Rev Ecol Evol Syst. 2010;41:127–47.
Kroeker K, Kordas R, Crim R, Singh G. Meta-analysis reveals negative yet variable effects of ocean acidification on marine organisms. Ecol Lett. 2010;13:1419–34.
Hofmann G, Washburn L. SBC LTER. Ocean: Time-series: Mid-water SeaFET and CO2 system chemistry at Mohawk Reef (MKO), ongoing since 2012-01-11. Santa Barbara Coastal LTER. 2015. doi:10.6073/pasta/826b170f29458104621aa9f0e36c8901.
Evans TG, Chan F, Menge BA, Hofmann GE. Transcriptomic responses to ocean acidification in larval sea urchins from a naturally variable pH environment. Mol Ecol. 2013;22:1609–25. doi:10.1111/mec.12188.
Wong JM, Johnson KM, Kelly MW, Hofmann GE. Transcriptomics reveal transgenerational effects in purple sea urchin embryos: Adult acclimation to upwelling conditions alters the response of their progeny to differential pCO₂ levels. Mol Ecol. 2018, 27(5):1120–1137. doi:10.1111/mec.14503.
Strader ME, Wong JM, Hofmann GE. Ocean acidification promotes broad transcriptomic responses in marine metazoans: a literature survey. Front Zool 2020, 17(7). doi:10.1186/s12983-020-0350-9.
Kültz D. Evolution of the cellular stress proteome: from monophyletic origin to ubiquitous function. J Exp Biol. 2003;206:3119–24.
Johnson K, Hofmann G. Transcriptomic response of the Antarctic pteropod Limacina helicina antarctica to ocean acidification. Bmc Genomics. 2017, 18:812. doi:10.1186/s12864-017-4161-0.
Moya A, Huisman L, Forêt S, Gattuso JP, Hayward DC, Ball EE, Miller DJ. Rapid acclimation of juvenile corals to CO2-mediated acidification by upregulation of heat shock protein and Bcl-2 genes. Mol Ecol. 2015;24(2):438–52. doi:10.1111/mec.13021.
Carreiro-Silva, Cerqueira T, Godinho A, Caetano M, Santos RS, Bettencourt R. Molecular mechanisms underlying the physiological responses of the cold-water coral Desmophyllum dianthus to ocean acidification. Coral Reefs. 2014;33:465–76. doi:10.1007/s00338-014-1129-2.
Lindquist S. The heat-shock response. Annu Rev Biochem. 1986;55:1151–91.
Sørensen JG, Kristensen TN, Loeschcke V. The evolutionary and ecological role of heat shock proteins. Ecol Lett. 2003;6:1025–37. doi:10.1046/j.1461-0248.2003.00528.x.
Wernberg T, Smale DA, Thomsen MS. A decade of climate change experiments on marine organisms: procedures, patterns and problems. Glob Chang Biol. 2012;18:1491–8. doi:10.1111/j.1365-2486.2012.02656.x.
Qiu J, Wang W-N, Wang L-j, Liu Y-F, Wang A-L. Oxidative stress, DNA damage and osmolality in the Pacific white shrimp, Litopenaeus vannamei exposed to acute low temperature stress. Comp Biochem Phys C. 2011;154:36–41.
Sorte CJB, Hofmann GE. Thermotolerance and heat-shock protein expression in Northeastern Pacific Nucella species with different biogeographical ranges. Mar Biol. 2005;146:985–93. doi:10.1007/s00227-004-1508-2.
Clegg JS, Uhlnger KR, Jackson SA, Cherr GN, Rifkin E, Friedman CS. Induced thermotolerance and the heat shock protein–70 family in the Pacific oyster Crassostrea gigas. Mol Mar Biol Biotechnol. 1998;7(1):21–30.
Brun NT, Bricelj VM, MacRae TH, Ross NW. Acquisition of thermotolerance in bay scallops, Argopecten irradians irradians, via differential induction of heat shock proteins. J Exp Mar Biol Ecol. 2009;371:77–83.
Dong Y, Dong S. Induced thermotolerance and expression of heat shock protein 70 in sea cucumber Apostichopus japonicus. Fish Sci. 2008;74:573–8.
González K, Gaitán-Espitia JD, Font A, Cárdenas CA, González-Aravena M. Expression pattern of heat shock proteins during acute thermal stress in the Antarctic sea urchin, Sterechinus neumayeri. Rev Chil Hist Nat 2016, 89. doi:10.1186/s40693-016-0052-z.
Hammond LM, Hofmann GE. Thermal tolerance of Strongylocentrotus purpuratus early life history stages: mortality, stress-induced gene expression and biogeographic patterns. Mar Biol. 2010;157:2677–87. doi:10.1007/s00227-010-1528-z.
Clark MS, Suckling CC, Cavallo A, Mackenzie CL, Thorne MAS, Davies AJ, Peck LS. Molecular mechanisms underpinning transgenerational plasticity in the green sea urchin Psammechinus miliaris. Sci Rep. 2019, 9.
Krebs RA, Feder ME. Hsp70 and larval thermotolerance in Drosophila melanogaster: how much is enough and when is more too much? J Insect Physiol. 1998;44:1091–101.
Evans TG, Watson-Wynn P. Effects of seawater acidification on gene expression: Resolving broader-scale trends in sea urchins. Biol Bull. 2014;226(3):237–54.
Padilla-Gamiño J, Kelly M, Evans T, Hofmann G. Temperature and CO₂ additively regulate physiology, morphology and genomic responses of larval sea urchins, Strongylocentrotus purpuratus. Proc R Soc B 2013, 280. doi:10.1098/rspb.2013.0155.
Stumpp M, Hu MY, Melzner F, Gutowska MA, Dorey N, Himmerkus N, Holtmann WC, Dupont ST, Thorndyke MC, Bleich M. Acidified seawater impacts sea urchin larvae pH regulatory systems relevant for calcification. PNAS. 2012;109(44):18192–7.
Todgham A, Hofmann G. Transcriptomic response of sea urchin larvae Strongylocentrotus purpuratus to CO₂-driven seawater acidification. J Exp Biol. 2009;212:2579–94. doi:10.1242/jeb.032540.
Martin S, Richier S, Pedrotti M-L, Dupont S, Castejon C, Gerakis Y, Kerros M-E, Oberhänsli F, Teyssié J-L, Jeffree R, et al. Early development and molecular plasticity in the Mediterranean sea urchin Paracentrotus lividus exposed to CO₂-driven acidification. J Exp Biol. 2011;214:1357–68. doi:10.1242/jeb.051169.
Kelly MW, Padilla-Gamiño JL, Hofmann GE. Natural variation and the capacity to adapt to ocean acidification in the keystone sea urchin Strongylocentrotus purpuratus. Glob Chang Biol. 2013;19:2536–46. doi:10.1111/gcb.12251.
Pespeni M, Chan F, Menge B, Palumbi S. Signs of adaptation to local pH conditions across an environmental mosaic in the California Current Ecosystem. Integr Comp Biol. 2013;53(5):857–70.
Pespeni M, Sanford E, Gaylord B, Hill T, Hosfelt J, Jaris H, LaVigne M, Lenz E, Russell A, Young M, et al. Evolutionary change during experimental ocean acidification. PNAS. 2013;110(17):6937–6924.
Kapsenberg L, Okamoto D, Dutton J, Hofmann G. Sensitivity of sea urchin fertilization to pH varies across a natural pH mosaic. Ecol Evol. 2017, 7:1737–1750. doi:10.1002/ece3.2776.
Evans TG, Pespeni MH, Hofmann GE, Palumbi SR, Sanford E. Transcriptomic responses to seawater acidification among sea urchin populations inhabiting a natural pH mosaic. Mol Ecol. 2017, 26(8):2257–2275. doi:10.1111/mec.14038.
Kroeker K, Sanford E, Rose J, Blanchette C, Chan F, Chavez F, Gaylord B, Helmuth B, Hill T, Hofmann G, et al. Interacting environmental mosaics drive geographic variation in mussel performance and predation vulnerability. Ecol Lett. 2016;19:771–9. doi:10.1111/ele.12613.
Hofmann G, Evans T, Kelly M, Padilla-Gamiño J, Blanchette C, Washburn L, Chan F, McManus M, Menge B, Gaylord B, et al. Exploring local adaptation and the ocean acidification seascape – studies in the California Current Large Marine Ecosystem. Biogeosciences. 2014;11:1–13.
Lester SE, Tobin ED, Behrens MD. Disease dynamics and the potential role of thermal stress in the sea urchin, Strongylocentrotus purpuratus. Can J Fish Aquat Sci. 2007;64(2):314–23. doi:10.1139/f07-010.
Frölicher TL, Laufkötter C. Emerging risks from marine heat waves. Nat Commun. 2018, 9(650). doi:10.1038/s41467-018-03163-6.
Ummenhofer CC, Meehl GA.. Soc B. Extreme weather and climate events with ecological relevance: a review. Philos T R. 2017;372:20160135. doi:10.1098/rstb.2016.0135.
Byrne M, Foo SA, Soars NA, Wolfe KDL, Nguyen HD, Hardy N, Dworjanyn SA. Ocean warming will mitigate the effects of acidification on calcifying sea urchin larvae (Heliocidaris tuberculata) from the Australian global warming hot spot. J Exp Mar Biol Ecol. 2013;448:250–7. doi:10.1016/j.jembe.2013.07.016.
Byrne M, Lamare M, Winter D, Dworjanyn SA, Uthicke S. The stunting effect of a high CO₂ ocean on calcification and development in sea urchin larvae, a synthesis from the tropics to the poles. Philos Trans R Soc Lond B Biol Sci 2013, 368(1627). doi:10.1098/rstb.2012.0439.
Sheppard Brennand H, Soars N, Dworjanyn S, Davis A, Byrne M. Impact of ocean warming and ocean acidification on larval development and calcification in the sea urchin Tripneustes gratilla. PLoS One 2010, 5(6). doi:10.1371/journal.pone.0011372.
Feely R, Sabine C, Hernandez-Ayon J, Ianson D, Hales B. Evidence for upwelling of corrosive “acidified” water onto the continental shelf. Science. 2008;320:1490–2.
Pennington JT, Chavez FP. Seasonal fluctuations of temperature, salinity, nitrate, chlorophyll and primary production at station H3/M1 over 1989–1996 in Monterey Bay, California. Deep Sea Res Part 2. 2000, 47:947–973.
Bennett J, Giese AC. The annual reproductive and nutritional cycles in two western sea urchins. Mar Biol Lab. 1995;109(2):226–37.
Kato S, Schroeter SC. Biology of the red sea urchin, Strongylocentrotus franciscanus, and its fishery in California. Marine Fisheries Review. 1985;47(3):1–20.
Strathmann MF: Reproduction and Development of Marine Invertebrates of the Northern Pacific Coast. USA. University of Washington Press. 1987.
Bakun A, Black BA, Bograd SJ, García-Reyes M, Miller AJ, Rykaczewski RR, Sydeman WJ. Anticipated effects of climate change on coastal upwelling ecosystems. Current Climate Change Reports. 2015;1(2):85–93. doi:10.1007/s40641-015-0008-4.
Sydeman WJ, García-Reyes M, Schoeman DS, Rykaczewski RR, Thompson SA, Black BA, Bograd SJ. Climate change and wind intensification in coastal upwelling ecosystems. Science. 2014;345(6192):77–80.
Snyder MA, Sloan LC, Diffenbaugh NS, Bell JL. Future climate change and upwelling in the California Current. Geophys Res Lett 2003, 30(15). doi:10.1029/2003GL017647.
Caputi N, Kangas M, Denham A, Feng M, Pearce A, Hetzel Y, Chandrapavan A. Management adaptation of invertebrate fisheries to an extreme marine heat wave event at a global warming hot spot. Ecol Evol. 2016;6(11):3583–93. doi:10.1002/ece3.2137.
Wilson JR, Lomonico S, Bradley D, Sievanen L, Dempsey T, Bell M, McAfee S, Costello C, Szuwalski C, McGonigal H, et al. Adaptive comanagement to achieve climate-ready fisheries. Conserv Lett. 2018;11:e12452. doi:10.1111/conl.12452.
IPCC. Summary for Policymakers. 2019.
Krueger F. Trim, Galore!. A wrapper tool around Cutadapt and FastQC to consistently apply quality and adapter trimming to FastQ files. Available at: . 2015.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012;9:357–9.
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, et al. vegan: Community Ecology Package. R package version 25 – 6. 2019.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289–300.

Download PDF

Journal Publication

published 07 Jan, 2021

Read the published version in BMC Genomics →

Editorial decision: Major revision
27 Jul, 2020
Review #2 received at journal
26 Jul, 2020
Reviewer #2 agreed at journal
13 Jul, 2020
Review #1 received at journal
12 Jul, 2020
Reviewers invited by journal
07 Jul, 2020
Reviewer #1 agreed at journal
07 Jul, 2020
Editor assigned by journal
02 Jul, 2020
Submission checks completed at journal
29 Jun, 2020
Editor invited by journal
29 Jun, 2020
First submitted to journal
22 Jun, 2020

You are reading this older preprint version

Read the latest preprint version →

Gene Expression Patterns of Red Sea Urchins (Mesocentrotus Franciscanus) Exposed to Different Combinations of Temperature and pCO2 During Early Development

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Results

Summary statistics and overview of RNA-seq

Developmental stage influenced transcriptomic patterns

Temperature treatment affected prism gene expression

Gene expression modules were correlated with treatments and physiology

Discussion

Transcriptomic patterns varied by developmental stage: General patterns

Temperature influenced gastrula embryos on a molecular level

Cellular stress response

Transmembrane transport and metabolism

Regulation of gene expression

Metabolism and ion transport

Cellular stress response

Temperature was the dominant factor at the prism stage

Similarities to the transcriptomic response at the gastrula stage

DNA damage and gene expression regulation at lower temperatures

Temperature did not influence HSP gene expression

Performance under current and future ocean conditions

Conclusions

Methods

Animal collection and culturing

Sample generation and sequencing

Data processing

DE Analysis and WGCNA

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1