Genomic Characterization of the Evolutionary Potential of the Sea Urchin Strongylocentrotus droebachiensis Facing Ocean Acidification.

Ocean acidification (OA) is increasing due to anthropogenic CO2 emissions and poses a threat to marine species and communities worldwide. To better project the effects of acidification on organisms' health and persistence, an understanding is needed of the 1) mechanisms underlying developmental and physiological tolerance and 2) potential populations have for rapid evolutionary adaptation. This is especially challenging in nonmodel species where targeted assays of metabolism and stress physiology may not be available or economical for large-scale assessments of genetic constraints. We used mRNA sequencing and a quantitative genetics breeding design to study mechanisms underlying genetic variability and tolerance to decreased seawater pH (-0.4 pH units) in larvae of the sea urchin Strongylocentrotus droebachiensis. We used a gene ontology-based approach to integrate expression profiles into indirect measures of cellular and biochemical traits underlying variation in larval performance (i.e., growth rates). Molecular responses to OA were complex, involving changes to several functions such as growth rates, cell division, metabolism, and immune activities. Surprisingly, the magnitude of pH effects on molecular traits tended to be small relative to variation attributable to segregating functional genetic variation in this species. We discuss how the application of transcriptomics and quantitative genetics approaches across diverse species can enrich our understanding of the biological impacts of climate change.

Introduction

The rapid pace of climate change threatens the persistence of species worldwide by exposing individuals to environmental conditions outside the range of recent historical precedents (IPCC 2014). In the marine realm, the combination of elevated temperatures and decreased pH (i.e., ocean acidification [OA]) poses a particularly acute challenge (Caldeira and Wickett 2003; Doney et al. 2009). Phenotypic responses to environmental factors may be a key determinant of a species’ prospects for survival, by buffering critical life processes via altered behavior, phenology, or stress response mechanisms (Williams et al. 2008). Evolutionary adaptation provides another path for species persistence and may be particularly advantageous over the long term (Davis and Shaw 2001; Davis et al. 2005). The prospects for evolutionary rescue depend on many factors, including demography, life history, and the amount of standing genetic variation in traits directly related to the tolerance of the environmental challenge (Lynch and Lande 1993; Bürger and Lynch 1995; Gomulkiewicz and Holt 1995). Even if a population harbors extensive genetic variation that influences traits necessary for adaptation, genetic or developmental trade-offs can place limitations on its ability to adapt (Lande and Arnold 1983). Although many studies have documented genetic variation that influences the tolerance of stressors related to climate change, evidence for recent adaptation—particularly adaptation that can keep up with the rapid pace of climate change—remains limited (Gienapp et al. 2008; Hendry et al. 2008; Kelly et al. 2013 Vander Wal et al. 2013).

Most evolutionary and ecological studies focus on accessible traits such as growth, survival, or phenology (e.g., Parmesan and Yohe 2003; Hendry et al. 2008; Gienapp et al. 2008), whereas variation in physiological, metabolic, and developmental traits is less commonly considered even though they are also key determinants of environmental responses (Wikelski and Cooke 2006; Chown and Gaston 2008; Somero 2010; Pan et al. 2015). Importantly, the link between these internal traits and external metrics of performance is often not obvious (Houle 2010). A comprehensive, unbiased approach is therefore needed to understand the genetic basis for plasticity and evolutionary responses to climate change.

A central challenge facing such efforts is that key physiological and developmental traits are often difficult, time consuming, or expensive to measure, making it impractical to survey a population for genetic variation affecting these traits. Genomic approaches, particularly ones that combine high-throughput molecular phenotyping with measures of high-level performance traits, provide a partial but highly practical solution. Transcriptome-wide gene expression assays are increasingly affordable and produce a rich source of data, both on variation in gene expression and on variation in cellular, physiological, and developmental processes that leave signatures in the covariation of functionally related genes (Alvarez et al. 2015). To list just a few examples, recent studies using transcriptomics in evolutionary ecology have identified gene networks underlying adaptive developmental plasticity in cichlid fish (Schneider et al. 2014), physiological mechanisms of cold and hypoxia tolerance in deer mice (Cheviron et al. 2014), and hormone signals regulating drought responses in Arabidopsis (Marais Des et al. 2012). Analytical approaches that leverage databases of gene function or genetic pathways (gene set enrichment analyses; Subramanian et al. 2005) are particularly promising: they reduce the high-dimensional data to a more accessible scale, highlight variation in informative molecular processes, and increase the power to detect subtle signals in noisy gene expression data. Here, we apply a transcriptome-based systems genetics approach (Civelek and Lusis 2013) to study the molecular response to changes in ocean pH relevant in the context of OA in larvae of the green sea urchin, Strongylocentrotus droebachiensis.

In marine systems, OA is one of the most immediate results of climate change and is recognized as a major threat to many species and communities (IPCC 2014; Breitburg et al. 2015; Gaylord et al. 2015). CO2 from the atmosphere dissolves in the ocean and reacts with water, lowering seawater pH and altering the carbonate balance of seawater (Doney et al. 2009). OA can affect a wide range of physiological processes in marine species, including acid–base balance, metabolism, development, deposition of skeleton, and growth rates (Dupont and Pörtner 2013). A number of studies in recent years have investigated the tolerance of sea urchins to decreased pH (reviewed in Dupont et al. 2010; Byrne 2011). In general, sea urchins have proven surprisingly robust to low pH (Dupont and Thorndyke 2013). However, large differences in tolerance are apparent even among closely related species, different life history stages can differ in their tolerance of low pH, and effects may increase when combined with other climate-related stressors such as increasing temperatures (Byrne 2011; Gianguzza et al. 2014). At the molecular level, larvae grown under low pH show increased expression of genes related to metabolic and ion regulation, and differential regulation of several genes involved in stress responses, apoptosis, calcification and skeletal formation (Todgham and Hofmann 2009; Martin et al. 2011; Stumpp et al. 2011a; Evans et al. 2013; Padilla-Gamiño et al. 2013; Pespeni et al. 2013). Several studies have also observed genetic variation in the response of growth rates (Sunday et al. 2011; Kelly et al. 2013) and development (Foo et al. 2012) to decreased pH in sea urchin embryos and larvae, as well as genetic changes in allele frequencies in artificial populations grown under simulated future conditions (Pespeni et al. 2013).

We extend these studies using transcriptome-wide gene expression profiling and a quantitative genetic breeding design to quantify molecular plasticity and test for segregating genetic variation affecting cellular and physiological traits in cultures of larvae spawned from S. droebachiensis individuals collected in northern Europe. We demonstrate that the molecular response to decreased pH in these larvae is broad and complex. Surprisingly, we also find that the gene expression plasticity to pH is dwarfed in magnitude and frequency by the influence of genetic variation among families. Finally, we show how combining signals over many functionally related genes can provide detailed information regarding 1) molecular mechanisms underlying responses to decreased pH and 2) aspects of these mechanisms that are mostly likely to contribute to near-term adaptations.

Materials and Methods

Experimental Design

Traits with ample genetic variation can efficiently respond to natural selection. A traditional tool to assess genetic variation in a population is to test for consistent differences among the offspring with different fathers (Lynch and Walsh 1998) because sperm is thought to contribute little to zygotes beyond DNA although the generality of this conclusion has been challenged by recent studies (Krawetz 2005; Curley et al. 2011; Crean et al. 2013). Consistent differences among offspring with different mothers may also be caused by genetic (heritable) variation, but may additionally be caused by maternal environmental variation that is passed transgenerationally through egg quality or epigenetic modifications. We used a three-way factorial breeding design to quantify effects of genetic background and seawater pH on larval gene expression, survival, and growth. The three factors were as follows: father (seven levels), mother (two levels), and pH (two levels; 8.0 and 7.6). A decrease by 0.4 units from the average pH at the study site (pHT 8.0) is within the range of the projections from the IPCC (2014) for average pH in 2100. Additionally, taking into account local environmental variability, a pH of 7.6 is the extreme of present environmental variability (Dorey et al. 2013). Cultures of larvae representing all 42 combinations of these three sets of factors were created and monitored for up to 9 days.

Animal Collection and Rearing

Adult S. droebachiensis were collected during the fall of 2009 from northern Norway (three males and one female; Hammerfest: 70°39’N/23°39’E, collected by SCAN-AQUA) and northern Denmark (four males and two females; Anholt: 56°42N/11°31’E, collected by the Marine Biological Laboratory of Helsingør). The original experimental design intended to test for differences between the populations in early-life responses to low pH; however, we do not focus on these differences here. Adults were maintained in a deep-water flow-through aquarium at the Sven Lovén Centre for Marine Sciences - Kristineberg (Fiskebäckskil, Sweden) and were fed ad libitum using Ulva spp. for 6–7 months before starting the experiment (March 2010). Spawning was induced by injecting ∼2 ml 0.5 mM KCl into the coelomic cavity. Sperm (four males from Denmark and three from Norway) was collected by pipette and kept dry until fertilization. Eggs (two females from Denmark and one from Norway) were collected in 250-ml beakers of filtered ambient seawater (FSW; 0.22 μm) and washed once. Total egg number was estimated by counting small-volume samples of the egg solution, and ∼200K–350K eggs were divided among seven 250-ml beakers to be fertilized, each with sperm from a different male. Fertilization was followed by monitoring the elevation of the vitelline envelope under a stereomicroscope (fertilization in all cultures was >95%). This created 21 families of full-sib embryos (fig. 1).

. 1.—

Experimental layout. We used a three-way factorial layout to assess the effects of pH (control vs low pH), male parent (seven individuals), and female parent (three individuals). Cultures were created representing all 42 combinations of these three factors. Larvae were raised for 14 days, and monitored daily for mortality, growth rate (days 1−8), and seawater chemistry (table 1).

Table 1

Seawater Carbonate Chemistry in Cultures at Each pH Treatment

Treatment	Control (n = 41)	Low pH (n = 96)
pHT	8.00 ± 0.03	7.62 ± 0.11
T (°C)	10.17 ± 0.20	10.34 ± 0.24
DIC (mmol/kg)	2.00 ± 0.08	2.13 ± 0.06
pCO2 (μatm)	419 ± 39	1145 ± 296
Ωar	1.76 ± 0.14	0.83 ± 0.25
Ωca	2.77 ± 0.22	1.30 ± 0.40

Note.—Seawater total scale pH (pHT), temperature (T), and DIC were measured daily in four control bottles and eight to 20 randomly-chosen acidified bottles (10 days). These measurements were used to calculate CO2 partial pressure (pCO2) as well as aragonite and calcite saturation states (respectively Ωar and Ωca), assuming a salinity of 34.7, using the package seacarb for R. All the values are expressed as mean ± SD.

Each family of embryos was split between two 5-l Erlenmeyer flasks of FSW for culturing at an initial density of 3–5 embryos/ml. All culture vessels were preequilibrated at the target pH (table 1. control: ca. pHT ≈ 8.00, ≈420 μatm, and low pH: pH T ≈ 7.6, ≈1,150 μatm). Culture vessels were aerated and gently mixed with a stream of bubbles of pressurized air and kept in the dark. Cultures were monitored for growth and survival until sampling.

Seawater pH, temperature, and dissolved inorganic carbon (DIC) were monitored throughout the experiment (table 1). Culture pH was maintained (±0.04 units) using Aqua Medic pH controllers (NBS scale, AquaMedic, Germany) that controlled valves for injection of pure CO2 into the culture flasks once pHNBS exceeded the target value. pH-system settings were adjusted from pH measurements on the total scale (pHT) using TRIS (Tris/HCl) and AMP (2-aminopyridine/HCl) buffer solutions with a salinity of 32.0 (provided by Unité d'Océanographie Chimique, Université de Liège, Belgium). DIC was measured with an automated carbon dioxide analyzer (CIBA Corning 965 UK). pCO2 and calcium carbonate saturation states for calcite and aragonite (Ωca and Ωar) were calculated from DIC and pHT using the R package seacarb (Gattuso et al. 2015).

Larval Mortality

On day 6, two 10-ml samples of each culture were aspirated and counted. This provided an estimate of the concentration of surviving larvae in each culture. Larval mortality associated with father or pH treatment was tested by comparing larval concentrations among the cultures with the same mother using a three-way analysis of variance (ANOVA) because the batch of eggs from each mother was distributed equally among all of her cultures. The ANOVA table is presented in supplementary table S1, Supplementary Material online.

Larval Growth Rate

Average growth rates were estimated for each culture over the first 7 days. Each day, 10 larvae were randomly sampled from each culture, fixed with 4% PFA, and imaged with a Leica microscope mounted with a DFC295 camera. From these images, the average body length (see measurement from Stumpp et al. 2011b) of each culture on each day was calculated. Growth rates were calculated in mm/day using linear regression of body length on age.

Gene Expression Measurements

Gene expression was measured using RNAseq on pools of young unfed embryos from each of the 42 cultures. Collection times were calibrated for each culture individually based on the per-culture growth rates and varied between 5.4 and 9.4 days when the average larval body length was approximately 0.36 mm. Sampling based on size rather than on physical time accounts for potential effects of pH on developmental rates (Pörtner et al. 2010; Stumpp et al. 2011a, 2011b). Developmental stages were consistent at the time of sampling, with most larvae in all cultures just starting to form buds for the posterodorsal pair of skeletal arms (supplementary figure S1, Supplementary Material online).

Approximately 3,000 larvae were collected by gently pouring 1–l of each culture into a 1.5-cm sieve consisting of a submerged nylon mesh and then pipetting the concentrated larvae into a clear 2-ml screw-top microcentrifuge tube. Once in the tube, the larvae were quickly mixed again and two 100-μl samples of highly concentrated larvae (50–200 larvae each) were collected and fixed in 4% PFA to estimate the number of collected larvae and the average larval size. Remaining larvae were pelleted by centrifugation, the seawater removed by aspiration, and 600 ml of RLT buffer (Qiagen, Hilden, Germany) added to stabilize the RNA. Tubes were vortexed rapidly for 15 s to lyse the cells and then stored at −80°C. The entire process from sieve to RLT took about 1 minute. RNA was extracted with the RNeasy kit (Qiagen) after treatment with DNase, and RNA quality was assessed with a total RNA analysis ng sensitivity (Eukaryote) assay on an Agilent 2100 Bioanalyzer. No samples showed signs of RNA degradation. The resulting mRNA was purified by poly-A selection.

RNAseq libraries were prepared by the Duke GCB Sequencing and Genomic Technologies Shared Resource using the Illumina TruSeq v1 kit with poly-A selection to purify mRNA. Six samples were pooled per lane to be run on the Illumina HiSeq 2000 system with version 3 chemistry. 1.3 × 109 50-bp single-end reads were generated, representing 1.8–5.3 × 107 reads per sample (median 3.2 × 107). Because no genome sequence of S. droebachiensis is available, we mapped reads to the genome sequence of the closely related species, Strongylocentrotus purpuratus (v3.1; Sodergren et al. 2006). We used the program bowtie2 (v2.0.0-beta5; Langmead et al. 2009; Langmead and Salzberg 2012) to map reads with the “–local” option and settings “-D 20 -R 3 -N 1 -L 20 -i S,1,0.50.” Across individuals, 74–86% of reads mapped to the S. purpuratus genome, and 67% of these mapped uniquely. Only uniquely mapped reads were used to estimate gene expression.

To convert mapped reads into gene expression measurements, we used the program HTseq-count (Anders et al. 2015) to count the number of reads that aligned uniquely to any exon, or either the 5' or 3’ untranslated region of any isoform of the 29,016 v3.1 GLEAN gene models (GLEAN-UTR-3.1.gff3, http://www.echinobase.org/Echinobase/SpDownloads; last accessed July 26, 2012). Gene expression analyses were performed in R (R Development Core Team 2010). Sample normalization factors were calculated with the calcnormfactors function of the edgeR package (Robinson and Oshlack 2010; Robinson et al. 2010) using default parameters. Genes with fewer than 10 counts over the 42 samples were not analyzed further, leaving a total of 22,430 genes.

Gene Set Variation Analysis

We used gene set variation analysis (GSVA; Hänzelmann et al. 2013) implemented in the R package GSVA to explore variation in higher-order molecular traits. GSVA defines a set of synthetic traits based on a list of gene ontology terms. Each trait’s value represents the extent to which genes labeled with a specific term tend to be up- or downregulated in that sample, measured using a Kolmogorov–Smirnov-like random walk statistic. Coordinated regulation of genes in the same pathway or involved in similar processes is evidence for changes in the activity of specific cellular functions (Subramanian et al. 2005).

We used two sources for gene function assignments: the PANTHER database (Mi et al. 2013) and a sea urchin-specific database of 130 ontology terms annotated by the sea urchin research community ("hand-annotated ontology": Sodergren et al. 2006; Tu et al. 2012). The PANTHER ontology is a reduced set terms from the Gene Ontology Project. For the sea urchin, the PANTHER ontology (PTHR8.1_sea_urchin) consists of 128 molecular function (MF) terms, 165 biological process (BP) terms, 22 cellular component (CC) terms, and 171 protein class (PC) terms with at least 10 and not >1,000 genes. We report results for all classes of ontology terms in the Supplementary Material online, but only discuss results from the MF and hand-annotated ontologies because we find them to be the most straight-forward to interpret. The hand-annotated ontologies are less comprehensive than the PANTHER ontologies, but include annotations for processes that have been especially studied in sea urchins, including biomineralization-related genes. Phenotypic correlations among gene set traits within and among ontologies can help guide interpretation and are presented as heatmaps in supplementary figure S1, Supplementary Material online.

We used the bubble plot method described by Supek et al. (2011) to visualize relationships among gene set traits. Within each of the PANTHER ontologies, specific terms are associated in a hierarchical tree, with more general terms (containing more genes) toward the base of the tree and more specific terms (containing fewer genes) toward the tips. We calculated the SimREL Semantic Similarity (Schlicker et al. 2006) distance between terms in the ontology based on the tree distance and term sizes and represent this distance on the 2D plane of the bubble plot using multidimensional scaling (R function sammon).

Statistical Analyses

All statistical analyses were performed in R v3.2.2 (R Development Core Team 2010).

Our experimental design consisted of three sets of factors: father, mother, and pH treatment fig. 1. The goal of the statistical analysis was to identify traits affected by each of these factors. For each trait, we considered a factor to be important if the main effect or any of the interactions of that factor with other factors were significant. Therefore, our lists of pH-responsive traits include all traits affected by pH in the offspring of any father or any mother. Traits are considered to be genetically variable if father effects were detected in either pH treatment or offspring of any of the three mothers.

To quantify the importance of each factor, we used the R Bioconductor package limma (Ritchie et al. 2015) to fit the following Gaussian linear models to the variation in each trait: where Y was the trait value of the sample in pH treatment i with father j and mother k, pH, F, and M are effects of the specific levels of each factor, PopF and PopM are effects of the source population (Norway or Denmark) of each parent, terms with x’s denote interactions, and μ is a global intercept. The source population terms accounted for any systematic differences (genetic or environmental) between urchins collected from each population. Although genetic differences among populations can serve as a source of variation for adaptation, trans-generational environmental effects cannot. To be conservative, we excluded differences in trait means or responses to pH between source populations from our estimates of the amount of genetic variation in each trait. Given our limited sample size (3–4 fathers and 1–2 mothers per population), we were unable to effectively test for population differences in any trait.

For the gene expression data, we converted counts to log2CPM (counts-per-million, with prior.count = 1) with the function cpm (Robinson et al. 2010), calculated sample-weights (accounting for variation in sample quality) with the arrayWeights function (Ritchie et al. 2006) and then fit model (1) and calculated smoothed Type III F-statistics for the sets of coefficients associated with each factor using the functions lmFit and eBayes (Smyth 2004) using the limma-trend method (Law et al. 2014). The sample-weights calculated for the gene expression data were carried over to the analysis of the gene set variation traits because these relied on the same gene expression values, but not to the analyses of survival and growth rate traits. Empirical Bayes smoothing of SEs was performed separately for the gene expression traits, each class of gene set traits, growth rate, and the survival traits.

We calculated the percentage of variance explained by a factor as the total weighted sum-of-squares for that factor (including interactions) after accounting for all other factors divided by the total weighted sum-of-squares of the observations around the global mean. We used the qvalue package (Storey et al. 2015) to calculate family-wise false discovery rates (FDR) and estimated numbers of true positives (1−π0) for each set of tests. We tested for enrichments of genes with functional annotations among genes with small P values using the Wilcoxon rank-sum test (Wilcox.test).

We used multiple linear regressions to test for associations between gene expression variation and variation in performance. For each performance trait (growth rate or survival fraction), we ran a penalized regression of the trait value on either the entire gene expression matrix (22,430 genes), or one of the classes of gene set traits using the LASSO penalty (Tibshirani 1996) implemented in the R function glmnet (Friedman et al. 2010). For each regression, we used leave-one-out cross-validation to select the tuning parameter lambda and scored models based on the improvement in the mean squared error of prediction.

Results

All trait values are presented in supplementary table S2, Supplementary Material online, and statistics for all tests of pH and parental effects are presented in supplementary table S3, Supplementary Material online.

Impact of Seawater pH

We did not observe any effects of low pH seawater on larval mortality (supplementary table S1, Supplementary Material online). Low pH seawater caused a 14% reduction in larval growth rate (low pH: 0.033 mm/day, control: 0.039 mm/day; P < 0.0001, fig. 2).

. 2.—

Larval growth rates are reduced in low pH seawater. Boxplots show median and quartiles of the distributions of growth rates for cultures grown in control or low pH seawater (N = 21 for each). The effect of pH treatment on growth rate was significant (P = 2.78 × 10−8). Growth rates (mm/day) were calculated using linear regression of daily measures of ∼10 larvae/culture over the first 6–9 days of development.

Gene expression responses to low pH seawater were common but generally subtle. Using the qvalue function (Storey et al. 2015), we estimated that 43.4% of the 22,430 measured genes were affected by the pH treatment in any of the family groups (cultures sharing a father or a mother), but the average response to low pH across all cultures was smaller than 2-fold for all but nine genes (log2FC < +/- 1, fig. 3, supplementary table S1, Supplementary Material online). Overall, among the 517 genes, we could declare significant at an FDR of 5% (corresponding to a P value threshold of 0.002, fig. 3, supplementary table S1, Supplementary Material online), many more genes were upregulated than downregulated at low pH (65%).

. 3.—

Physiological and molecular responses attributed to seawater pH are common but subtle relative to variation attributed to father or mother effects. (A) Mean (across all cultures) log2 responses to low pH seawater for 22,430 genes plotted against mean log2(expression). Genes with significant responses to low pH seawater (5% FDR) are highlighted in blue. (B) Boxplots show median and quartiles of the distributions of percentage of total among-culture variation in each gene expression trait (n = 22,430) or each MF gene set (n = 128) trait that could be attributed to each experimental factor (mother, father or pH).

In order to gain insight into biological responses to low pH, we used gene ontology annotations to group genes into functionally similar modules and measured how these modules responded to the pH treatment. The set of gene ontology annotations we used included annotations for 14,951 genes, which were significantly enriched for responses to the pH treatment (Wilcoxon sign-rank W: 5.49 × 107, P = 0.016). Using GSVA (Hänzelmann et al. 2013), we combined the expression levels of these 14,951 genes into measures of 616 synthetic gene set traits, divided into five categories (PANTHER ontologies: MF: 128 terms, BP: 165 terms, CC: 22 terms, PC: 171 terms, and hand-annotated ontology: 130 terms). We focused on MF and hand-annotated categories that are the easiest to interpret. Results for all gene set traits are presented in supplementary table S1, Supplementary Material online. pH responses in the 128 MF gene set traits were larger in magnitude (relative to the total among-sample variance, Wilcoxon sign-rank W: 1.73 × 106, P < 0.0001) and were more common than the responses to pH of the individual genes (fig. 3). We estimated that 60.9% of these MF gene set traits were affected by the pH treatment. pH treatment effects (either as main effects or as interactions with genetic backgrounds) accounted for >25% of the among-sample variation in 38 of these traits and >50% of the variation in 2. The bubble plot in figure 4 summarizes the 25 MF gene set traits that we could declare to be significantly affected by low pH in any cohort at a 5% FDR. Several key groups of MF gene sets constitute the core of the response to pH. Downregulated gene sets were mostly related to DNA and RNA metabolism and cell division. Upregulated gene sets were related to cell signaling and other transmembrane functions such as ion channel activities and G-protein coupled receptor activity.

. 4.—

Summary of the response to pH in Molecular Function gene set traits. Bubble plot of MF gene set variation traits with significant (5% FDR) responses to low pH seawater. Bubble plots represent: (i) The a priori relationship among the MF gene sets. Bubble centers are arranged based on a multidimensional scaling projection (R function sammon) of the SimREL distances (Schlicker et al. 2006) among the PANTHER MF ontology terms. This distance takes into account both the tree-relationships among terms and the number of genes in each category. (ii) The number of genes linked to each MF term (bubble area is proportional to gene number). (iii) The percentage of variation in each gene set explained by the pH treatment and the mean direction of the response to low pH (orange = increase in expression, blue = decrease in expression). This plot is based on the REVIGO scatter plot (Supek et al. 2011).

Results were largely congruent for the 130 gene set traits based on the hand-annotated ontologies. Low pH seawater affected 71.2% of these traits and explained >25% of the variance in 46. The 62 hand-annotated gene set traits we could declare to be significantly affected by low pH in any cohort primarily involved processes such as the cell cycle, cell–cell signaling and immune-related processes (supplementary table S3, Supplementary Material online). The category of biominerlization genes was significantly downregulated by low pH (−0.204).

Genetic and Maternal Effects on Larval Traits

Across both performance and molecular traits, differences associated with larval parentage tended to be larger and more common than differences induced by low seawater pH (fig. 3). We detected significant mother effects on growth rate (P = 0.001), but no significant father effects (P = 0.44). Among the 22,430 genes, 77.4% had expression variation associated with inherited father effects and 48.9% had expression variation associated with mother effects (note that the experiment included seven fathers and only three mothers, but only among parent differences within the same source population were counted). Again, genes with functional annotations were enriched for significant father or mother effects relative to unannotated genes. Gene expression profiles integrated into gene set traits also showed high levels of parent-associated variation: 93.8% of the 128 MF gene set traits showed father-associated variation and 55.7% showed mother-associated variation.

Interactions between Genetic Variation and the Response to Low pH

By splitting cultures of genetically related larvae between the two pH treatments, we were able to test for differences in responses to pH among culture with different fathers or mothers, and differences in the father effects when crossed to different mothers. We estimated that 20% of genes had a different expression response to pH across families with different fathers and 6.7% of genes had a different response across cultures with different mothers. Among the 128 MF gene set traits, 1% had a different response to pH among cultures with different mothers, but we found no evidence for a different response to pH in any of these traits among cultures with different fathers. We estimated that 6.7% of genes and 28.6% of MF gene set traits had different father effects depending on the mother. However, due to the low level of replication of these interactions in our experimental design, we were not able to declare more than a handful to be significant while maintaining a FDR < 5% (supplementary table S3, Supplementary Material online) and the true frequencies of these gene–environment, or gene–gene interactions may be higher.

Relationship between Molecular Traits and Larval Performance

To test if the molecular traits could help explain variation in larval performance, we used multiple regression with a LASSO penalty to try to predict the variation in growth rate among cultures based on the whole matrix of 22,430 genes, or each of the five classes of gene set traits. We used the LASSO penalty because it provides variable selection among the predictors and helps with model regularization when there are more predictors than samples. The BP and hand-annotated gene set traits were the best predictors of growth rate variation, as measured by cross-validation performance, and out-performed the raw gene expression data (table 2). The best model selected 25 BP gene set traits and improved the mean squared error by 79.2% relative to a model with only culture pH. These 25 BP gene set traits are displayed in figure 5. In this model, the BP traits pyrimidine base metabolic process (positive) and RNA localization (negative) were selected as the most strongly associated with growth rate.

Table 2

Percent Decrease in Mean Squared Error of Prediction (MSE) Values Relative to pH Treatment Alone for the Regression of Growth Rate on Each Set of Molecular Traits

Molecular Trait Class	Growth Rate
BP	79.2% (24)
MF	46.8% (13)
CC	18.5% (9)
PC	36.9% (29)
Hand	64.7% (41)
Gene expression	56% (40)

Note.—MSE was calculated using the cv.glmnet function with alpha = 1 for the LASSO penalty. Values represent 1−MSE(full)/MSE(pH) for each model at the optimal value of the lambda tuning parameter. The number in parentheses is the number of molecular traits with non-zero regression coefficients in the best model.

. 5.—

Biological process gene set traits associated with variation in larval growth rates. Bubble plot show the 22 BP gene set traits selected by the LASSO penalized regression as important predictors of growth rate variation. Color hue and intensity represents the sign and magnitude of each of the non-zero regression coefficients in the best model. The lambda parameter was chosen by leave-one-out cross validation.

Discussion

Molecular Responses to Low pH

Our results demonstrate that the molecular response to low pH in larvae of S. droebachiensis is wide-ranging and complex. We observed plasticity in response to low pH in the expression of nearly half of all assayed genes, and in three quarters of the higher-order molecular traits inferred by integrating the expression variation across groups of functionally related genes. Previous studies investigating the effects of low pH on gene expression in sea urchin larvae have reported changes in the expression of genes related to calcification and biomineralization (Todgham and Hofmann 2009; Martin et al. 2011; Stumpp et al. 2011a; Padilla-Gamiño et al. 2013; Evans et al. 2013), cytoskeleton and cell division (Todgham and Hofmann 2009; Padilla-Gamiño et al. 2013), acid/base regulation (Todgham and Hofmann 2009; Stumpp et al. 2011a; Stumpp et al. 2015), and metabolism (Todgham and Hofmann 2009; Stumpp et al. 2011a), with the direction and strengths of the responses varying by species and the particulars of the environment treatments (e.g., simultaneous increases in temperature). We observed changes in many of these same processes in S. droebachiensis larvae. This is consistent with physiological studies showing that S. droebachiensis pluteus larvae are unable to compensate for an extracellular acidosis (pHe) resulting from an exposure to low pH. However, the calcifying primary mesenchyme cells are able to fully compensate an induced intracellular acidosis (pHi) using a bicarbonate buffer mechanism involving secondary active Na+-dependent membrane transport proteins (Stumpp et al. 2012). Additional energetic costs also derived from compensatory mechanisms associated with larval gastric pH changes (Stumpp et al. 2013). These extra costs lead to a shift in energy budget, with less energy available for growth and leading to a delay in development (Dupont and Thorndyke 2013; Jager et al. 2016). Our results also highlighted a strong signal of increased expression of many genes involved in immune responses and other responses to extracellular environments, perhaps indicating an involvement of the immune system in pH tolerance or an increased sensitivity of the larvae to infections or additional environmental stresses. This is in contrast with studies of mussels (Bibby et al. 2008) and starfish (Hernroth et al. 2011) both showed a depressed immune system under low pH but consistent with an increased in cellular immune response in two adult echinoderm species, including S. droebachiensis (Dupont and Thorndyke 2012). Dupont and Thorndyke (2012) hypothesized that there may be a direct link between pHe and cellular immune-response. As it was demonstrated that decreased environmental pH also lead to an uncompensated pHe decrease in pluteus larvae, the observed upregulation of the immune system may result from such a link. However, other measures of cellular stress did not show evidence of a pH response: no consistent expression changes were observed in the class of heat shock proteins, or other molecular chaperones, unlike what has been seen in other environmental stressors, such as temperature (Hammond and Hofmann 2010; Runcie et al. 2012).

The increased number of significant effects in our data relative to previous studies is largely a function of our increased sample size (21 cultures per treatment), as the magnitude of pH effects on gene expression were generally small as previously reported (Todgham and Hofmann 2009; Pespeni et al. 2013). Despite this molecular plasticity, we found a general robustness of larval development to predicted near-future levels of pH, with the larvae appearing healthy despite slower growth rates and no increase in larval mortality due to pH. This replicates similar findings in previous studies (Dupont and Thorndyke 2008; Martin et al. 2011; Stumpp et al. 2011b). In the absence of data from long-term cultures, however, potential impacts of lower pH on survival to metamorphosis and lifetime fecundity remain unclear.

Natural Variation in Larval Traits

Despite clear evidence of physiological, developmental, and growth impacts of pH on the larvae, the pH treatments did not account for a particularly large portion of the total gene expression variation among cultures of larvae. pH treatment effects accounted for <25% of the among-culture variation in the majority of the gene expression traits measured (fig. 3). Instead, our results show that genetic variants and/or effects of differences in paternal or maternal environments contributed to a much larger fraction of the trait variation (fig. 3). Although it is not possible using our experimental design to unequivocally differentiate between paternal environmental effects and genetic variation, a high frequency and magnitude of genetic variation in gene expression has been observed before in sea urchins (Runcie et al. 2012; Garfield et al. 2013), as well as in other taxa including fish, plants, Drosophila, and primates (Oleksiak et al. 2002; Schadt et al. 2003; Nuzhdin et al. 2008; Ayroles et al. 2009). We extended these results by demonstrating that there is a high level of genetic variation in higher-order molecular traits, captured by the gene set traits defined by ontology terms from the PANTHER gene ontology (Mi et al. 2013), and by the sea urchin-specific gene function ontology (Sodergren et al. 2006; Tu et al. 2012). For example, the father and mother effects on the MF signatures were as frequent and of as large a magnitude (proportionally to the total amount of trait variation) as the parental effects on the expression of individual genes.

Prospects for Evolutionary Rescue

The large amount of genetic variation in molecular gene expression and survival fractions suggests that rapid adaptation to low pH seawater might be possible in this species. Large amounts of standing genetic variation can accelerate adaptation if the variation contributes to variation in fitness, and if genetic correlations among traits with contrasting effects on fitness are not too strong (Walsh and Blows 2009). We observed considerably more genetic variation (father-associated variation) in the molecular traits than in larval growth rates (fig. 3). It is possible that genetic variation can accumulate in low-level traits like gene expression because gene networks are structured to buffer much of this variation from affecting high-order phenotypes such as morphology (Runcie et al. 2012; Garfield et al. 2013). However, at least a portion of the variation in the gene expression traits was correlated with variation in growth rate (fig. 5 and table 2) and thus may be relevant for fitness. It is also possible that this molecular variation may be a source of "cryptic" genetic variation that could be exposed by changing environments (Gibson and Dworkin 2004), and thus accelerate adaptation in changing environments. We found evidence for gene and environment interactions being important for a small percentage of the gene expression traits, but our sample sizes were too small to clearly identify any patterns as to which traits might show a change in genetic variance in acidified seawater. The fact that we observed as strong evidence for genetic variation in the higher order gene set traits as in the individual genes suggests that much of this genetic variation may be highly pleiotropic. The possibility of extensive pleiotropy is an important caveat in the interpretation of our finding that many genes show variation in expression. This is because pleiotropy places genetic constraints on the direction of evolutionary change (Walsh and Blows 2009), and thus the capacity of a population to adapt in any particular way may be more limited than the magnitude of genetic variation in each trait initially suggests.

Another limitation of our finding is the use of the two scenarios used in this study: the average pH at present (pHT 8.0) and the average pH projected for 2100 (pHT 7.6). However, both pH values are within the present range of natural variability (Dorey et al. 2013). Previous physiological works on the impact of acidification across a wide range of pH on the development of S. droebachiensis have demonstrated a clearly pH-dependent increase in stress with decreasing pH, including an apparent tipping point around a pH of 7.5 (Jager et al. 2016). Below that threshold, different physiological processes were involved leading to increased mortality at low pH (Dorey et al. 2013). Different evolutionary processes may also be involved, as demonstrated in copepods. When exposed to decreased pH with the present day’s range of natural variability, copepods responded mostly through phenotypic plasticity. When exposed to pH below present day’s natural variability, transgenerational effects (including genetic adaptation) set in (Thor and Dupont 2015; De Wit et al. 2015).

Gene Expression Profiling in Climate Change Research

The use of RNAseq to address problems in ecological and evolutionary functional genomics is growing rapidly and shows great potential to greatly improve our understanding of the key traits underlying adaptation and acclimation to changing climates (Gracey 2007; Whitehead 2012; Alvarez et al. 2015; Papetti et al. 2016; Todd et al. 2016). Transcriptome profiling by RNAseq can now be applied to nearly any organism. The principal advantage of RNAseq for climate change research is that it rapidly profiles a wide range of molecular processes, many of which are difficult to measure directly (Evans and Hofmann 2012). For example, techniques to measure respiration and immune system activity are technically challenging to apply to large numbers of individuals (Hernroth et al. 2011; Stumpp et al. 2011b). With gene expression profiling, we can quickly sort through hundreds of conceivable molecular processes to identify ones that are likely relevant and are worth further study. As sequencing costs continue to decline rapidly, it is becoming feasible to measure gene expression in many hundreds or even thousands of samples. Such sample sizes are necessary for accurately estimating genetic variances and genetic correlations among traits that are necessary for predicting evolutionary change (Roff 1995).

By themselves, transcriptome-wide gene expression data can be overwhelming in their scale. Individual gene expression measurements are often too noisy (whether due to measurement error or systematic variation such as genetic effects) to use as reliable signatures of underlying molecular trait variation. Directly incorporating information on gene function (Hänzelmann et al. 2013) and modeling the correlation among traits (Runcie and Mukherjee 2013) can dramatically increase both the power to detect molecular variation, and make the results easier to interpret. There are important limitations: characterizing traits based on functional classification of genes is inherently imprecise would benefit from explicit information about how genes in gene sets interact (Khatri et al. 2012) and requires follow-up studies to directly measure the traits identified by gene expression signatures. Nevertheless, by simultaneously measuring many aspects of physiology and development, gene expression can provide a vastly more comprehensive estimate of the potential for adaptation to climate change (De Wit et al. 2015; Rose et al. 2016).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.