Little evidence of adaptation potential to ocean acidification in sea urchins living in “Future Ocean” conditions at a CO2 vent

Abstract Ocean acidification (OA) can be detrimental to calcifying marine organisms, with stunting of invertebrate larval development one of the most consistent responses. Effects are usually measured by short‐term, within‐generation exposure, an approach that does not consider the potential for adaptation. We examined the genetic response to OA of larvae of the tropical sea urchin Echinometra sp. C. raised on coral reefs that were either influenced by CO2 vents (pH ~ 7.9, future OA condition) or nonvent control reefs (pH 8.2). We assembled a high quality de novo transcriptome of Echinometra embryos (8 hr) and pluteus larvae (48 hr) and identified 68,056 SNPs. We tested for outlier SNPs and functional enrichment in embryos and larvae raised from adults from the control or vent sites. Generally, highest F ST values in embryos were observed between sites (intrinsic adaptation, most representative of the gene pool in the spawned populations). This comparison also had the highest number of outlier loci (40). In the other comparisons, classical adaptation (comparing larvae with adults from the control transplanted to either the control or vent conditions) and reverse adaptation (larvae from the vent site returned to the vent or explanted at the control), we only observed modest numbers of outlier SNPs (6–19) and only enrichment in two functional pathways. Most of the outliers detected were silent substitutions without adaptive potential. We conclude that there is little evidence of realized adaptation potential during early development, while some potential (albeit relatively low) exists in the intrinsic gene pool after more than one generation of exposure.

Most tests of OA on marine invertebrates, however, have been conducted on individuals not acclimatized to these altered experimental conditions, and understanding the potential for adaptation of organisms is a major knowledge gap Riebesell & Gattuso, 2015). In fact, there are a number of studies that point to the potential for adaptation to reduced pH across a range of invertebrate taxa, for example, in bivalves (Zhao et al., 2018), crustaceans (Thor & Dupont, 2015) and echinoderms (Kelly, Padilla-Gamiño, & Hofmann, 2013). Similarly, experimentation on taxa living in naturally occurring CO 2 vent communities can increase our understanding of the acclimation potential of species both within  and across generations .
Next-generation sequencing of individuals or populations allows novel insights into population structures and the potential for lineages to adapt to potential environmental stressors over the longer term.
Tools for the analysis of RNA or DNA samples from pooled populations makes these tasks more time and cost effective (Kofler, Pandey, & Schlötterer, 2011;Schlötterer, Tobler, Kofler, & Nolte, 2014). In organisms with short generation times (e.g., Drosophila, Tobler et al., 2013), "Evolve and Resequencing" studies can be conducted, where changes in allele frequencies for a large number of single nucleotide polymorphisms (SNPs) and outlier analysis can identify which genes are under adaptive pressure, and how many generations would be required for adaptation . Similarly, physiological responses can be measured after exposing many generations to altered OA conditions (Sunday et al., 2014). In the marine realm, this can be achieved for short-generation taxa such as coccolithophorids, where growth under different climate change conditions was observed over 460 generations (Schluter et al., 2014).
In marine studies to date, pool-seq approaches have been mainly applied in fishes. For instance, stickleback (Gasterosteus aculeatus) populations in the Baltic Sea exhibited genetic differentiation which was consistent with local adaptation to temperature and salinity gradients (Guo, DeFaveri, Sotelo, Nair, & Merilä, 2015). A similar approach has been applied to Herring populations (Clupea harengus) which exhibited many outlier loci with apparent adaptation to local environmental conditions (Guo, Li, & Merilä, 2016). RNA-seq approaches can also be applied to gene expression studies, for example, to investigate changes in heat stress tolerance in Copepods (Schoville, Barreto, Moy, Wolff, & Burton, 2012).
For populations of the sea urchin Strongylocentrotus purpuratus, Pespeni, Garfield, Manier, and Palumbi (2012)  The study identified significant allelic changes in more than 40 functional protein classes and the authors concluded that this sea urchin has a high capacity for fast adaptation to OA conditions and that "…standing genetic variation could be a reservoir of resilience to climate change." Follow-up studies also used a pool-seq approach for SNP-based population genomics (Pespeni, Chan, et al., 2013;Pespeni, Sanford, et al., 2013) and differential expression studies (Evans, Chan, Menge, & Hofmann, 2013;Evans, Pespeni, Hofmann, Palumbi, & Sanford, 2017) to investigate the potential of that species for adaptation to ocean acidification along environmental pH gradients and in controlled laboratory experiments.
Sea urchins (Echinodermata: Echinoidea) are considered "keystone species" in many marine ecosystems and have important ecological functions as habitat modifiers, especially through macroalgal grazing (Lessios, 2001). For example, die-off of a Caribbean species has resulted in overgrowth of coral areas with algae and an ecosystem phase-shift (Hughes, 1994). In contrast, increased urchin abundances resulting from predator overfishing can lead to increased bioerosion and to reduced coral reef recovery after disturbances (Carreiro-Silva & McClanahan, 2001). Increased urchin abundances can also cause "deforestation" of temperate kelp forest, ultimately resulting in urchin barrens (Steneck et al., 2002).
Experimentation on the effects of ocean acidification on sea urchins have been relatively extensive, with the response in planktonic larvae more distinct among species and ecosystems compared with those seen in adults . The tropical coral inhabiting sea urchin genus Echinometra is relatively well-studied. The effects of OA on adult Echinometra growth depend on species and experimental setup. A strong decline in growth with increasing Pco 2 was reported for a Caribbean species (Courtney, Westfield, & Ries, 2013), and Echinometra species often demonstrated a reduction in growth of the larval skeleton with increased Pco 2 and more asymmetric larvae (Kurihara & Shirayama, 2004;Lamare et al., 2016;Uthicke, Soars, Foo, & Byrne, 2013). By contrast, the effects on several Indo-Pacific species were only subtle (Uthicke, Liddy, Nguyen, & Byrne, 2014;Uthicke, Soars, et al., 2013) or were not detected (Hazan, Wangensteen, & Fine, 2014;Moulin et al., 2015; Grosjean, Leblud, Batigny, & Dubois, 2014).
The effects of OA on developmental stages of marine invertebrates, such as sea urchins, are usually investigated in relatively shortterm (weeks to months) controlled aquarium experiments. A recent alternative to laboratory experiments, and which allows long-term exposures and observations, is to study organisms or ecosystems at CO 2 vents in the field. This also considers secondary ecological effects of OA (e.g., increased food availability due to enhanced algal or seagrass growth). Vents used for these studies now include Mediterranean rocky shores (Hall-Spencer et al., 2008), White Island New Zealand (Nagelkerken, Russell, Gillanders, & Connell, 2015), and Mariana Islands and Papua New Guinea (PNG) for coral reef environments (Enochs et al., 2015;Fabricius et al., 2011). We previously investigated sea urchin adult growth  and larval development  in Echinometra sp. C at the well-studied vent sites on Papua New Guinean coral reefs. Research in CO 2 vents can be exploited to explore the phenotypic and genotypic responses in Echinometra offspring of acclimated animals.
These vents release nearly 100% pure CO 2 gas into a coral reef systems and are seen as ideal "natural laboratories" to predict future responses to OA conditions (Fabricius et al., 2011). Surprisingly, for adult urchins we found increased growth rates at the vent sites. This was despite a low pH that is similar to that expected at the end of this century . By contrast, in situ Echinometra larval growth studies in these sites showed that exposure to future OA conditions can lead to stunted and abnormal growth . This is in agreement with previous laboratory studies on the same genus (Uthicke, Soars, et al., 2013) and on a wide variety of sea urchins from tropical to arctic regions .
Here, we exposed larvae in situ to conditions resembling OA and to present-day seawater (control) conditions at a CO 2 vent system in Papua New Guinea. We present a de novo assembled transcriptome of early embryos (8 hr) and 48 hr old pluteus larvae of Echinometra sp. C. Subsequently, we investigated the potential of adaptation of these urchins at the vent sites using a population genomic approach.
Using >60,000 SNPs and outlier analysis, we investigated genes and functional groups under apparent selection. We compared larvae derived from adult urchins from either a control or vent location that were cross transplanted to the reciprocal site with those exposed at their parent's origin ( Figure 1). This design should allow us to distinguish outlier SNPs in genes providing potential for adaptation to OA during short incubations, from those already differentiated in the two source populations.

| Species and study sites
Experiments were conducted on Echinometra sp. C sensu McCartney, Keller, and Lessios (2000) as evidence by cytochrome oxidase subunit 1 (CO1) sequencing . In situ larval experiments were conducted at two sites in Milne Bay Province, PNG (Fabricius et al., 2011). These consist of a CO 2 vent site (9.82426°S, 150.81788°E), and an adjacent ambient control site ~500 m from the vent (9.82807°S, 150.82027°E) at Upa-Upasina, NW Normanby Island (see Figure 1 in Lamare et al., 2016). The pH and Pco 2 at the control site (pH T : 7.92-7.92, Pco 2 : 410-450 ppm) corresponds to present-day conditions, whereas average pH and Pco 2 values at the vent site (pH T : 7.57-7.64 , , Pco 2 : 800-1,200 ppm) corresponds to future ocean acidification scenarios for the year ≈2100 under the RCP8.5 scenario (Meinshausen et al., 2011). There was no difference in water temperature during the deployments; average values for the control were 28.51, and for the vent 28.44°C , while previous studies did not detect a long-term difference in temperature between the control and vent site (Uthicke, Momigliano, et al., 2013). Echinometra sp. C collected from either Control (C) or Vent (V) sites, spawned, and their offspring explanted to both control locations (c) or vent locations (v). This gave rise to larvae with four histories (i.e., Cc = larvae from control adults explanted at control sites; Cv = larvae from control adults explanted at vent sites; Vc = larvae from vent adults explanted at control sites; Vv = larvae from vent adults explanted at vent sites). (b) Three potential Adaptive scenarios to changes in seawater pH for Echinometra sp. C larvae can be explored by detecting outlier SNPs between larvae of different histories. Classical adaptation (Cc vs. Cv) would indicate a potential for transgenerational adaptation at the genes carrying these outliers; Reverse adaptation (Vc vs. Vv) would indicate the potential for adaptation from ocean acidification to present-day conditions, and Intrinsic adaptation (Cc vs. Vv, representative of standing genetic variation) represent long-term (one full generation) intrinsic adaptation of the population at the vents. Classical adaptation and Reverse adaptation were tested after 8 and 48 hr of exposure. Intrinsic adaptation was only tested for 8 h larvae; SNPs in these larvae represent the closest approximation to allele frequencies in the parent populations from the control and the vent sites experimental run, adults (~25) were collected at either the control or vent site and spawned by injecting 0.5-1.0 ml of 0.5 M KCl. This procedure yielded 10 females and 8-9 males per experimental run.
Oocytes were cleaned and pooled, keeping the number of eggs per female approximately the same. Similarly, 1 µl of dry sperm from each male was mixed and the total activated in 10 ml filtered seawater and oocytes subsequently fertilized at a final concentration of 10 5 sperm/ ml. Following fertilization, embryos (700-900 individuals per replicate) for each experimental run were out-transplanted in replicated (N = 6) meshed 50 ml cages either to their original location or to a reciprocal location (i.e., embryos from control populations to vents sites, and vice versa). The larvae were then left to develop in situ for 48 hr, recollected, concentrated and immediately fixed in 8 subsamples of 2 ml RNAlater, stored at 4°C for 24 hr, then kept at −20°C until required. The 48h larvae were in the 2-4 arm pluteus stage.
To observe genetic responses in earlier life stages (i.e., up to gastrula), we prepared eight, 200 ml plastic jars with both control and vent water and added 700 to 900 embryos. These embryo cultures were incubated for 8 hr in 60 L flow-through aquaria to keep the temperature at ambient. At the end of the incubation, embryos were concentrated and fixed in 2 ml RNAlater as outlined above.

| RNA extraction and illumina sequencing
A subsample of each RNAlater fixed sample was counted, which indicated that all samples contained 200-500 larvae. Subsamples were combined and vacuum filtered onto a Durapore ® membrane filter (0.45 µm; Merck Millipore). Each filter was incubated in 1.5 ml TRIzol ® (Invitrogen) over night. After vortexing and centrifugation (4°C, 90 s, 10,000 g), 1.3 ml of the resulting supernatant was pipetted into a new tube and 260 µl chloroform added. Samples were then mixed, re-centrifuged (4°C, 15 min, 12,000 g) and 600 µl of the resulting supernatant transferred into a new tube. Samples were then carefully inverted after adding 600 µl of 70% ethanol. RNA extracts were subsequently purified using an RNeasy Mini Kit (Qiagen).
Purified RNA quality was checked by running all extracts on 2% agarose gels and determination of concentration and A260/A280 ratios on a Nanodrop™ 2000 (Thermo Fisher).
Libraries for sequencing were prepared using the TruSeq Stranded mRNA-seq kit (Illumina), using 1 µg of RNA as input for the poly-A pulldown followed by cDNA synthesis, A-tailing, and adapter ligation. The libraries were enriched by 15 PCR cycles and 75 bp paired-end sequencing was then done on an Illumina NextSeq 500 platform using the High Output v2 kit. The sequencing was carried out at the Ramaciotti Centre for Genomics at the University of New South Wales, Australia.

| De novo transcriptome assembly and assessment
The SolexaQA toolkit (http://solex aqa.sourc eforge.net/) was used to quality trim the raw reads. Quality trimmed reads from individual samples were assembled separately with the Trinity assembly program version r2012-06-08 (Grabherr et al., 2011) using default parameters. A Reference Annotation Based Transcript (RABT) assembly approach was adopted, where transcripts with acceptable orthology to the model sea urchin S. purpuratus proteins were retained. This restricted the population differentiation analysis to variations identified in highly annotated sea urchin-specific transcripts.
The assembled transcripts from individual assembles were then combined. CD-HIT (Li & Godzik, 2006), a program for clustering nucleotide and protein sequences, was used to identify the representative transcripts to be included in the composite representative larval transcriptome. A Trinity module was used to identify near full length/ full-length transcripts when compared to the S. purpuratus protein set as reference. The BUSCO v3 tool (Waterhouse et al., 2017) was used for assessment of transcriptome completeness using the representation of core-conserved genes from the eukaryotic datasets.
OrfPredictor (Min, Butler, Storms, & Tsang, 2005) was used to predict putative open reading frames from the assembled transcriptome.

| Mapping
Quality-filtered reads from each site were aligned to the assembled transcriptome separately using BWA 0.6.1 ) with a maximum edit distance of two for a maximum insert size of 500 for a read pair, and one as the maximum number of alignment for reads paired properly. The mapping results in SAM format were converted into BAM format using SAMtools 0.1.18 ) and filtered for a minimum mapping quality of 20. Duplicate reads were removed using Picard Tools (http://broad insti tute.github.io/picar d/). BAM files were then converted into "mpileup" format using SAMtools 0.1.18 with maximum of 1,000 reads at a position per BAM file.
To verify the identity of the target species, we downloaded the mitochondrial NADH dehydrogenase genes ND1 and ND2 of 16 Echinometridae species from the NCBI database. We then also extracted the sequence of these genes from our transcriptome assembly and used Mega10 to align sequences and to construct individual phylogenetic trees by the neighbor-joining method, without distance corrections. We used 1,000 individual runs to test bootstrap support for individual clades.

| Detection of selection footprints (outlier SNPs)
We followed best practice guidelines for Pool-seq analysis (Schlötterer et al., 2014) for variant detection analysis of pooled populations. The duplicate datasets for each larval sample were combined for the SNP detection analysis. We employed the following steps when using PoPoolation2 (Kofler et al., 2011) along with specific parameters. The original sync file was subsampled to a target coverage of 50 reads. SNPs, which are possibly differentiated as a result of selection, were then identified using PoPoolation2 with a minimum minor allele count of 20 and coverage between 50 and 200 within each sample. Two methods were used for outlier detection. First, for each pairwise comparison F ST values for individual SNPs were calculated using PoPoolation2. SNPs falling into the 99.5 percentile of the empirical distribution of each pairwise F ST were identified as potentially differentiated loci following an empirical outlier detection approach (Akey et al., 2010;Guo et al., 2016). Fisher's exact test was used to determine whether any of the differences in allele frequencies were statistically significant. The negative log p-values obtained here were used downstream to identify possible enrichments using the tool ErmineJ (Lee, Braynen, Keshav, & Pavlidis, 2005). Secondly, all the SNPs identified using PoPoolation2 were converted to genotype data; this was done by simulation in the PoPoolation2 environment. The output from this was exported in GenePop format, and the transcript-specific files were merged into a single GenePop file and subsequently converted into BayeScan format. The BayeScan tool (Foll & Gaggiotti, 2008), which aims at identifying candidate loci under natural selection from genetic data and uses differences in allele frequencies between populations, was employed to estimate the posterior probability that a given locus is affected by selection. The following parameters were used for the BayeScan run: Prior odds of 10 were used to identify the top candidates of selected loci, with a thinning interval of 10, following 20 pilot runs of 5,000 iterations each, and a burn-in length of 50,000. A higher level of stringency was introduced into the outlier detection process wherein only the SNPs identified as outliers by both methods were considered as truly differentiated loci. This should ensure a smaller possibility of false positives.
In a separate analysis, genes differentially expressed (DE) across specific conditions were identified using the tool DESeq2 (Love, Huber, & Anders, 2014) with FDR ≤ 0.05 as a threshold. Outlier SNP-containing genes identified above, which were also found to be differentially expressed, were removed to avoid possible influence of gene expression on our results (Pespeni, Sanford, et al., 2013).
This overlap was <10% for each comparison (Tables S1 and S2). A full analysis of the DE genes was not presented here because it was not the aim of the study and thus the sampling design chosen not ideal for differential expression analysis; however, all DE genes identified in the main comparisons are presented in Tables S3 and S4.

| Enrichment analysis
Gene ontology annotations were used to check if the outlier SNPs identified with high allelic changes show enrichments for specific biological functions. We implemented the gene score resampling method from the tool ErmineJ (Gillis, Mistry, & Pavlidis, 2010) to test for correlation between functional enrichment and changes in allele frequency. Mean score was used for multiple SNPs encountered in a gene. Statistical significance was determined by using 10,000 permutations, and corrected p-values for false discovery rate were computed by the Benjamini-Hochberg approach. In addition, we tested for enrichment in biomineralization genes as proposed by Pespeni, Sanford, et al. (2013) and Evans et al. (2013). We used a list of 348 biomineralization (BM) genes from Evans and Watson-Wynn (2014) to identify BM gene orthologs in our larval transcriptome. We then tested if outlier SNPs identified were over-represented in each of the five main comparisons using Fisher's exact test.

| Larval development
Echinometra sp. C larvae analyzed here were from a subset of experiments reported by Lamare et al. (2016) that were undertaken at control and vent sites measuring morphometric changes at elevated Pco 2 concentrations. To put the current experiment into context, we briefly summarize these previous findings here. Compared to controls, 48 hr larvae exposed at vents exhibited a subtle but significant reduction in size, and larvae had increased abnormal features, measured as arm-length asymmetry. Larvae produced from adults acclimatized to ocean acidification conditions (i.e., comparing Cv and Vv larvae as defined here) did not show greater resilience to these future ocean conditions, suggesting no detectable acclimation to reduced pH in offspring from vent adults in terms of early development.  Table 1.

| Reference Annotation Based Transcript (RABT) and de novo transcriptome assembly
The tool BUSCO was used to explore the completeness of the transcriptome. It reported the presence of 277 (91%) Benchmarking Universal Single-Copy Orthologs (BUSCO) out of the 303 eukaryotic orthologs used in the assessment. Out of these, 190 orthologs were found to be in single copy, 87 in multi copies, 12 as fragments, and 14 were missing. This indicates that the assembled transcriptome is of high quality.
Clustering of de novo assembled ND1 and ND2 genes against reference sequences from GenBank confirmed that the assembled ND1 and ND2 genes clustered with other "Echinometra" species for both genes ( Figure S2) and most closely with another specimen of Echinometra sp. C.  Table 2). All overlap outliers identified were "high" outliers, thus indicative of diversifying selection (Figure 4). We also tested whether the substitutions in the outlier SNPs would lead to amino acid changes and thus provide potential for adaptation (Table 3, details in Table S3).

| Genome-wide population differentiation
Only 5 out of 66 (7.6%) outliers in the 8 hr data and 2 out of 21 (9.5%) of the 48 hr outliers led to amino acid changes, the remainder being silent substitutions. Several of these genes were uncharacterized, and the remaining 5 genes are not clustered in any biochemical pathway.
We tested for functional enrichment for individual sets of genes containing outlier SNPs for each comparison for biological pathways using ErmineJ. We did not detect any enrichment for all but one comparison in the 48h adaptation scenario ( Table 2). In that comparison, we found two pathways enriched. These pathways were similar, namely "regulation of protein localization in plasma membrane" (GO:1903076) and "positive regulation of protein localization to plasma membrane" (GO:1903078). We identified 274 genes in the larval assembly as putative orthologs of genes involved in biomineralization, however, in only two of our five main comparisons did an outlier SNP occur in a biomineralization gene (Table 2). Hence, none of the comparisons had a significant enrichment of outliers in these genes.

| D ISCUSS I ON
This research tested the hypothesis that long-term exposure of adult sea urchins to elevated Pco 2 (intrinsic adaptation) or short-term (8, 48 hr) larval exposure (classic adaptation) can lead to adaptation to ocean acidification conditions. We will argue that we must reject this hypothesis because: (a) compared to other studies the number of outliers detected is small; (b) most of these outliers are "silent" substitutions, and; (c) apparent outliers exhibit little enrichment for biological pathways.
Given that Echinometra sp. C is a nonmodel species and the only urchin genome available is from a different echinoid Family (Strongylocentrotus, Sodergren et al., 2006), we assembled a transcriptome for Echinometra sp. C from all available larvae samples. A high N50 suggested that the transcriptome was highly contiguous in its structure. The completeness of the transcriptome was demonstrated by using BUSCO, which showed that 91% of core-conserved genes are found in the assembly. The number of unique proteins identified in our study was comparable to recently assembled sea ur-  (reviewed in: Byrne et al., 2013), including larvae of Echinometra spp. (Uthicke, Soars, et al., 2013). However, to our knowledge none of these studies investigated mortality; thus, it is possible that measured stress in these early phases does not translate into (selective) mortality. The stress could nevertheless affect fitness and capacity to grow into adulthood, and however, this was not investigated here.   (Uthicke & Benzie, 2003), sea stars (Harrison, Pratchett, Messmer, Saenz-Agudelo, & Berumen, 2017;Williams & Benzie, 1997), or sea urchins (Duchaud, Durieux, Coupé, Pasqualini, & Ternengo, 2019;Tourón et al., 2018).
Also, SNP-based studies generally show low population genetic structure (F ST often < 0.002) within regions for planktotrophic echinoderms such as holothurians (Xuereb et al., 2018). Based on mitochondrial genes, intra-regional differences between Echinometra sp. C can be nonsignificant (bindin gene) or range from 0.04 to 0.15 (COI, Geyer & Palumbi, 2003). An allozyme study on Echinometra matthaei demonstrated F ST values around 0.01 on geographic scales from single to 1,000's of kilometers (Watts, Johnson, & Black, 1990); thus, F ST values detected here over small distances may not be unusual but also reflect genetic patchiness.
Outliers SNPs in the intrinsic comparison were 2-4 times more common than in the other comparisons, indicating that some selection between the sites is realized. The outlier number detected in this comparison was in the same order as in an experimental study on larvae of S. purpuratus (~40, Pespeni, Sanford, et al., 2013) but much lower than found in adults of the same species along natural pH gradient (>300, Pespeni, Chan, et al., 2013).
A reduction of size and increased abnormality upon exposure to ocean acidification is documented in many arctic to tropical sea urchin species . Morphological studies from the same experiments presented here  showed that, compared to larvae exposed at controls (Cc or Vc), 48h pluteus larvae exposed at vents (Cv or Vv) exhibited a reduction in size, and larvae had increased abnormality. Offspring from adults from the vents (i.e., those acclimatized to ocean acidification; V) were equally affected by OA when larvae were back-transplanted to the vents (Vv) compared to larvae from adults at the controls exposed at the vents (Cv). Thus, in terms of development, no evidence for selection or transgenerational phenotypic plasticity was detected in that study. Similarly, in an aquarium experiment on Echinometra sp. A (Uthicke, Soars, et al., 2013), larvae produced from adults pre-acclimatized (7 weeks) to low pH were not more resilient to the effects of ocean acidification compared to larvae with parents exposed to present-day conditions. Evans et al. (2017) found genes in tricarboxylic acid cycle, electron transport chain, and beta oxidation pathways upregulated in adult sea urchin larvae exposed to low pH (or high Pco 2 ). The authors concluded that this was an indication of ATP channeled toward maintaining acid-base homeostasis to increase tolerance to low pH, most likely at the cost of other measures of fitness. Given adults from the vent populations studied here grew faster and had a similar amount of energy channeled toward gonad production as those at controls, we saw no evidence for this in Echinometra sp. C adults . However, given Echinometra sp. larvae are generally smaller under OA and have more abnormalities Uthicke, Soars, et al., 2013) it is possible that these changes also occur, and without these, the effects of OA would be more drastic or lethal at earlier stages.
We discovered very few outliers indicative of potential ad-  (Evans et al., 2017).
None of the genes with nonsilent substitutions were differentially expressed on day 7 of Evan's study or were among the genes putatively under selection in other ocean acidification studies on S. purpuratus (Pespeni, Chan, et al., 2013;Pespeni, Sanford, et al., 2013).
Low statistical power to detect outlier loci in our study, especially for the intrinsic comparisons between two populations already separated by high F ST values may be one reason for detecting less evidence for selection then Pespeni, Sanford, et al. (2013). However, we used nearly twice as many parents as Pespeni et al., and the number of larvae sequenced (>200) was well in the range recommended for "Evolve and Resequencing" studies (Schlötterer et al., 2014).
Thus, at least for the "classical" and "reverse" comparison (using the same parents each, i.e., F ST at the start of the experiment close to 0), power cannot be assumed a reason for different findings compared to Pespeni et al. Even in the "intrinsic" comparison, we detected several apparent outliers with low minor allele frequencies (<0.1, Tables S1-S4), indication that low power may also not be a major issue in that comparison.
Based on previous studies on a temperate sea urchin, it was surprising that outlier genes in our study were not enriched for the biomineralization pathway. After 7 days of exposure to OA conditions, larvae of S. purpuratus outliers were significantly enriched in the biomineralization pathway (Pespeni, Sanford, et al., 2013), and this pathway was also over represented in outlier SNPs along a natural pH gradient (Pespeni, Chan, et al., 2013). However, there was no enrichment of differentially expression genes in the biomineralization pathway for larvae of the former species after 1 and 7 days exposure to OA (Evans et al., 2017).
Mortality rates may depend on the duration of experiments. This may explain the lower potential for adaptation observed in the present study on Echinomtera, compared to S. purpuratus (Pespeni, Chan, et al., 2013;Pespeni, Sanford, et al., 2013). In this respect, while both studies examined SNP outliers in the same developmental stage (i.e., 4-Armed pluteus), the timing of sampling was shorter in the present study (48 hr) compared with sampling time in (Pespeni, Chan, et al., 2013;Pespeni, Sanford, et al., 2013). It is possible, therefore, that the longer sampling duration may have allowed greater selection for genotype and hence greater differences in allele frequency seen in S. purpuratus. Although a previous study exhibited high (~50%) mortality in sea urchin larvae after 48 hr under OA (Gianguzza et al., 2014) it is possible that in our study, this length of time was insufficient to drive differential mortality or a resulting difference in allele frequency.
Furthermore, a potential for adaptation also requires the intrinsic presence of genetic variation leading to gene variants more suitable for altered environmental conditions. Some S. purpuratus populations inhabiting the west coast of North America naturally experience low pH in upwelling areas (Pespeni, Chan, et al., 2013), and possibly over many generations. Thus, standing genetic variance in that species may allow adaptation/acclimation to low pH on short time scales.
Although Echinometra sp. on shallow coral reef flats can be exposed to low pH (<7.8) and high Pco 2 (>1,000 ppm) at night (Uthicke et al., 2014), these short-term exposure to future ocean conditions may not be a sufficient evolutionary force resulting in potential for rapid adaptation and selection across generation.
In conclusion, our study revealed little evidence of a potential of adaptation during the first 2 days of larval development. Although some adaptation may have occurred in the intrinsic gene pool (separate populations at vents and controls), a previous morphometric study suggested that lifetime acclimation of the parents conveyed no resilience to OA to the larvae . Thus, the results from the present genetic study and the previous morphological study are congruent. The transcriptome assembly presented here will be an important resource for future adaptation studies of this keystone coral reef species.

ACK N OWLED G M ENTS
We thank the owners of the reefs in Upa-Upasina and Dobu for allowing us access to their reefs and to the councilors of the

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
SU has designed the study, ran experiments, lead the writing and analysis. ND has performed most of the data analysis and contributed to writing. MLi ran field experiments and contributed to writing. FP conducted laboratory analysis and contributed to writing.
MLa and MRW contributed significantly to the writing and data analysis writing.

DATA AVA I L A B I L I T Y S TAT E M E N T
The cleaned short read sequences and the transcriptome have been deposited in the European Nucleotide Archive (ENA