Cryptic genetic variation underpins rapid adaptation to ocean acidification

Global climate change has intensified the need to assess the capacity for natural populations to adapt to abrupt shifts in the environment. Reductions in seawater pH constitute a conspicuous stressor associated with increasing atmospheric carbon dioxide that is affecting ecosystems throughout the world’s oceans. Here, we quantify the phenotypic and genetic modifications associated with rapid adaptation to reduced seawater pH in the marine mussel, Mytilus galloprovincialis. We reared a genetically diverse larval population in ambient and extreme low pH conditions (pHT 8.1 and 7.4) and tracked changes in the larval size and allele frequency distributions through settlement. Additionally, we separated larvae by size to link a fitness-related trait to its underlying genetic background in each treatment. Both phenotypic and genetic data show that M. galloprovincialis can evolve in response to a decrease in seawater pH. This process is polygenic and characterized by genotype-environment interactions, suggesting the role of cryptic genetic variation in adaptation to future climate change. Holistically, this work provides insight into the processes underpinning rapid evolution, and demonstrates the importance of maintaining standing variation within natural populations to bolster species’ adaptive capacity as global change progresses.


Introduction 38
A fundamental focus of ecological and evolutionary biology is determining if and how 39 natural populations can adapt to rapid changes in the environment. Recent efforts that have 40 combined natural population censuses with genome-wide sequencing techniques have shown that 41 phenotypic changes due to abrupt environmental shifts oftentimes occur concomitantly to 42 signatures of selection at loci throughout the genome 1-4 . These studies demonstrate the 43 4 extent to which current levels of genetic variation will facilitate the magnitude and rate of 66 adaptation necessary for species persistence is unclear 20,21 . 67 In marine systems one pertinent threat is ocean acidification, the global-scale decline in 68 seawater pH driven by oceanic sequestration of anthropogenic carbon dioxide emissions 22 . The 69 current rate of pH decline is unprecedented in the past 55 million years 23 , and lab-based studies 70 have shown negative effects of expected pH conditions on a range of fitness-related traits (e.g., 71 growth, reproduction, and survival) across life-history stages and taxa 24 . Marine bivalves are one 72 of the most vulnerable taxa to ocean acidification 25,26 , particularly during larval development 27 . 73 The ecologically and economically valuable Mediterranean mussel, Mytilus galloprovincialis, is 74 an exemplary species for studying the effects of ocean acidification on larval development. Low 75 pH conditions reduce shell size and induce various, likely lethal, forms of abnormal larval 76 development 28,29 . Sensitivity to low pH, however, can vary substantially across larvae from 77 distinct parental crosses, suggesting that standing genetic variation could fuel an adaptive 78 response to ocean acidification 29 . 79 Here, we explored the potential for, and dynamics of, rapid adaptation to ocean 80 acidification in M. galloprovincialis. We tracked the phenotypic distributions and trajectories of 81 29,400 single nucleotide polymorphisms (SNPs), from the embryo stage through larval pelagic 82 growth and settlement in a genetically diverse larval population, within artificially imposed 83 ambient (pH T 8.05) and extreme low pH treatment (pH T 7.4). To test for a genotype-environment 84 interaction underpinning adaptation, we associated shell size, a trait negatively correlated with 85 fitness-related abnormalities 29 , to its underlying genetic background in each pH treatment. The 86 results presented demonstrate the capacity for natural populations to adapt to rapid 87 environmental change, and suggest that this process will be fueled by cryptic genetic variation. 88

Methods Summary 89
We quantified the effects of low pH exposure on phenotypic and genetic variation 90 throughout development in a single population of M. galloprovincialis larvae (Fig. 1). Larvae 91 were reared in ambient and low pH and (i) shell size distributions were quantified on days 3, 6, 92 7, 14, and 26; (ii) SNP frequencies across the species' exome were estimated on days 6, 26, and Replicate buckets marked with an "X" were destructively sampled (i.e. all larvae removed/preserved) and thus absent from the experimental system on subsequent sampling days. *Allele frequency data from two replicate buckets in the ambient treatment was generated on day 43, as the third replicate bucket to optimize protocol for sampling settled larvae  Results 111

Phenotypic Trajectories 112
As expected, shell size was significantly affected by pH treatment throughout the 113 experiment (linear mixed effects model: p = 0.029), and shell length of low pH larvae was 8% 114 smaller than that of larvae reared in ambient pH on days 3 and 7. Shell length was affected by the 115 interaction of day and treatment (linear mixed effects model: p < 0.001), indicating treatment 116 specific growth patterns. From days 7 to 26 the size distributions in each treatment began to 117 converge, with larvae in low pH becoming only 2.5% smaller than those cultured in the ambient 118 treatment by day 26 (Fig. 2). 119 120 121 122

Changes in Genetic Variation 123
We identified 29,400 SNPs across the species exome that were present within the larval 124 population across all sampling days and treatments. To link the observed phenotypic trends in 125 Figure 2 | Larval size distributions throughout the shell growing period. Larval size was significantly affected by treatment (p = 0.029) and the interaction of day and treatment (p < 0.001) throughout the shell growing period. stage and settlement (excluding size-separated groups) (Fig. 3). At later days (e.g. days 26 and 139 43), there is an observed increase in Euclidian distance among samples. This may be driven, at 140 least in part, by selection-induced declines in larval population survival throughout the pelagic 141 phase, and an associated increase in the influence of allele frequency "drift" among replicate 142 buckets. Observations of increased larval mortality throughout the experiment (indicated via 143 empty D-veliger shells in buckets) corroborated these trends, though we were unable to quantify 144 larval mortality (MCB, pers. obs.). 145 We identified SNPs that changed significantly in frequency (i.e., outlier SNPs) between 146 the day 0 larval population and the larvae sampled on day 6, 26, and 43 in each treatment. 147 Outlier SNPs were identified using a rank-based approach and the observed allele frequency shift 148  significantly changing in frequency throughout the course of the experiment relative to day 0 151 (Fig 4a). As the larvae sampled on day 6 were drawn from different replicate buckets as those 152 sampled on days 26 and 43 (Fig. 1), outlier SNPs observed on all three sampling days point to 153 candidate loci that may be putatively under selection in each pH environment. We used these 154 robust SNPs to generate lists of pH-specific loci or overlapping loci (genes containing outlier 155 SNPs that were responsive in each treatment). In total, we identified 99 ambient pH-specific loci 156 The size-separation of larvae on day 6 isolated the largest 18% from the smallest 82% in 180 the ambient and the largest 21% from the smallest 79% in the low pH treatment (Fig. S1). 181 Hereafter, these groups will be referred to as the fastest and slowest growers, respectively. Shell 182 size on day 6 was significantly affected by treatment and size class (Linear Mixed Effects Model, 183 p < 0.001). PCA using allele frequency data from the day 0 starting larval population and larval 184 samples collected on day 6 revealed a strong genetic signature of size class (Fig. 5). Specifically, 185 the fastest growers segregated along PC1 from the slowest growers in both treatments, with the 186 day 0 larval population and day 6 larval population samples (from each treatment) falling in 187 between the size-separated groups. The number of outlier SNP's differentiating the fastest and 188 slowest growers in each treatment, hereafter referred to as size-selected SNPs, was comparable: 189 Allele frequency data using 29,400 SNPs identified in larval smples collected on day 6, 6, and size separated larvae from day 6 (i.e. fastest and slowest growers).

Concurrent shifts in larval size and genetic variation throughout development 215
While previous work has shown strong negative effects of low pH on larval development 216 in bivalves 25,27 , the results presented here suggest that standing variation within the species could 217 facilitate rapid adaptation to ocean acidification. Observed shell length differences on days 3 and 218 7 matched expectations for the species based on previous work (-1 µm per 0.1 unit decrease in 219 pH) 29 . However, this difference was reduced ~50% by day 14. Mechanistically, low pH 220 treatment effects on bivalve larval shell growth are driven by the limited capacity of larvae to 221 regulate carbonate chemistry, specifically aragonite saturation state, in their calcifying space 31-33 . 222 Our data show that this physiological limitation is greatest prior to day 7, after which low pH 223 larvae were able to partially recover in size compared to larvae reared in the ambient treatment. 224 It is likely that partial recovery of shell size in low pH larvae observed by day 14 was, at 225 least in part, driven by natural selection. We have previously shown that the smallest D-veligers 226 in the low pH treatment display an increased prevalence of morphological abnormalities, which 227 likely become lethal during the shell growth period 29 . Directional selection against this 228 respectively. This suggests that a largely singular, intense selection event occurred prior to day 6 241 and may be responsible for the majority of genetic differentiation that occurs during larval 242 growth and settlement. We recently identified two specific early developmental processes that 243 are sensitive to low pH conditions and occur in this timeframe 29 . These processes include the 244 formation of the shell field (early trochophore stage) and the transition between growth of first 245 and second larval shell (late trochophore stage), both of which occur within 48 hours of 246 fertilization, resulting in a suite of size-dependent morphological abnormalities that likely 247 become lethal during the shell growth period 29 . Traditionally, metamorphosis from the 248 swimming D-veliger to the settled juvenile is regarded as the main genetic bottleneck during the 249 development of marine bivalve larvae 34 . Our sampling from the embryo stage through 250 13 settlement, however, suggests that there is a major selection event prior to day 6 that may have 251 an even larger effect on shaping genotypes of settled juveniles any selection thereafter. 252 Additional factors that may have allowed the larvae reared in led to the observed 253 phenotypic dynamics are food-augmented acclimation and selective mortality via food 254 competition. It has been demonstrated that increased energy availability can allow marine 255 invertebrates to withstand pH stress 35  The role of cryptic genetic variation in extreme low pH adaptation 272 We identified hundreds of loci responding to each pH treatment throughout the larval 273 period. Notably, 58% of our candidate low pH loci were statistically unchanged in the ambient 274 conditions. While some of this treatment disparity may be an artifact (i.e. false positives in the 275 low pH or false negatives in the ambient treatment), it is unlikely that this is the case for all the 276 unique SNPs identified. These data thus provide evidence that loci that were not critical for 277 fitness in the ambient environment came under strong selection in the low pH treatment. 278 Associating shell growth to its underlying genetic background in each environment strengthened 279 this conclusion. Specifically, there was an exceedingly small amount of overlap in size-selected 280 loci between treatments (76% of size-selected loci were unique to the low pH environment) and 281 a relative elevation of F ST differentiation between the fastest growers from the low and ambient 282 pH buckets. These patterns display a classic genotype-environment interaction in which a 283 particular genetic background exhibits a specific trait value in one environment (e.g. accelerated 284 shell growth in ambient pH), while an alternate genetic background leads to the same trait value 285 when the environment shifts 14 . As shell growth is a direct proxy for fitness 29,43 , these data 286 suggest that the most fit genotypes in ambient conditions may not be the individuals that harbor 287 the adaptive genetic variation necessary to improvs fitness in simulated ocean acidification. 288 Ultimately, our observation of genotype-environment interactions in response to a 289 dramatic shift in seawater pH provides strong evidence that adaptation to ocean acidification 290 may be fueled by cryptic genetic variation. The role of cryptic genetic variation during 291 adaptation to novel and extreme environmental conditions is becoming increasingly 292 recognized 12,14,15,44,45 , though its relative importance has yet to be demonstrated as clearly in 293 nature as it has been in theory 9,14,16 . Our data not only suggest the role of cryptic genetic 294 variation in rapid adaptation, but also demonstrate this phenomenon in the context of a non-295 model species subject to global change. The economic and ecological importance of marine 296 mussels, as well as their global exposure to declining seawater pH, highlights the need to 297 conserve standing variation in order to allow the adaptive capacity of natural populations to play 298 out as climate change progresses. 299 Many outlier loci in low pH were also outliers in the ambient treatment (42%). This

Putative targets of selection as ocean acidification progresses 313
The low-pH specific loci we identified (loci with outlier SNPs in every replicate bucket 314 and across all sampling days) provide targets of natural selection as ocean acidification 315 progresses. These included a HSPA1A gene (Swissprot ID: Q8K0U4), which encodes heat shock 316 protein 70 (HSP70), one of a group of gene products whose expression is induced by 317 physiological stressors and generally work to mediate/prevent protein denaturation and folding 48 318 . While substantial evidence has documented the role of HSP70 in the thermal stress response 319 across a range of taxa 48 , emerging transcriptomic studies have also demonstrated the protein's 320 role in the physiological response to low pH conditions in echinoderms 49 and bivalves 50 . Another 321 notable candidate locus is a tyrosinase-like protein (Swissprot ID: H2A0L0). Tyrosinase genes 322 are known to influence biomineralization in marine bivalves during larval 51,52 and juvenile 323 stages 53 , a process that is fundamentally affected by changes in seawater chemistry. While these 324 gene expression-based studies provided initial insight into the underlying physiological 325 responses to changes in seawater chemistry in marine bivalves, our study demonstrates the 326 presence of underlying genetic variation within these putative loci. This provides, to our 327 knowledge, the first documentation of standing genetic variation at functionally relevant loci 328 within marine bivalves, and ultimately offers robust evidence for the species' capacity to adapt to 329 extreme changes in seawater pH. We are currently investigating these candidates more intensely 330 through a combination of comparative transcriptomics, quantitative PCR, and in situ 331 hybridization in M. galloprovincialis larvae (Kapsenberg et al., unpublished). 332

Conclusions 333
Species persistence as global climate change progresses will, in part, hinge upon their 334 ability to evolve in response to the shifting abiotic environment 21 . Our data suggest that the 335 economically and ecologically valuable marine mussel, M. galloprovincialis, currently has the 336 standing variation necessary to adapt to ocean acidification, though with a potential trade-off of 337 shell size. We have demonstrated that much of this variation may be cryptic in ambient pH 338 conditions, yet bolsters a fitness-related trait (shell growth) when seawater pH is reduced. 339 Ultimately, these findings support conservation efforts aimed at maintaining variation within 340 natural populations to increase species resilience to future ocean conditions. In a broader 341 evolutionary framework, the substantial levels of genetic variation present in natural populations 342 have puzzled evolutionary biologists for decades 54 . Though this study does not address the 343 processes that maintain this variation during periods of environmental stasis, we demonstrate its 344 utility in rapid adaptation, thereby advancing our understanding of the mechanisms by which 345 natural populations evolve to abrupt changes in the environment. Within 3 weeks of the adult mussel collection, individuals were cleaned of all epibiota 356 using a metal brush, byssal threads were cut, and mussels were warmed in seawater heated to 357 27 o C (~+12 o C of holding conditions) to induce spawning. Individuals that began showing signs 358 of spawning were immediately isolated, and allowed to spawn in discrete vessels, which were 359 periodically rinsed to remove any potential gamete contamination. Gametes were examined for 360 viability and stored on ice (sperm) or at 16 o C (eggs). In total, gametes from 12 females and 16 361 males were isolated to generate a genetically diverse starting larval population. To produce 362 pairwise crosses, 150,000 eggs from each female were placed into sixteen separate vessels, 363 corresponding to the sixteen founding males. Sperm from each male was then used to fertilize 364 the eggs in the corresponding vessel, thus eliminating the potential effects of sperm competition 365 and ensuring that every male fertilized each female's eggs. After at least 90% of the eggs had 366 progressed to a 4-cell stage, equal volumes from each vessel were pooled to generate the day 0 367 larval population (~2 million individuals), from which the replicate culture buckets were seeded. 368 100,000 individuals were added to each culture buckets (N = 12, 18 embryos mL -1 ). The

Larval sampling 377
We strategically sampled larvae throughout the experiment to observe phenotypic and 378 genetic dynamics across key developmental events, including the trochophore to D-veliger 379 transition (day 6), the shell growth period (day 26), and the metamorphosis from D-veligers to 380 settlement (day 43). On day 6 of the experiment, larvae were sampled from three of the six 381 replicate buckets per treatment. A subset of larvae (N = 91-172) from each bucket was isolated to 382 obtain shell length distributions of larvae reared in the two treatments. The remaining larvae 383 (~25,000) were separated by shell size using a series of six Nitex mesh filters (70 µm, 65 µm, 60 384 µm, 55 µm, 50 µm, and 20 µm; Figure S1) and frozen at -80 o C. The smallest size group 385 contained larvae arrested at the trochophore stage, and therefore unlikely to survive. The 386 remaining five size classes isolated D-veligers from the smallest to the largest size. The shell 387 length distribution of the larvae was used to inform, a posteriori, which combination of size 388 classes would produce groups of the top 20% and bottom 80% of shell sizes from each treatment. 389 The relevant size groups from two replicates per treatment were then pooled for downstream 390 DNA analysis of each phenotypic group. For the third replicate, a posteriori, all size groups were 391 pooled in order to compute the allele frequency distribution from the entire larval population in 392 each treatment on day 6. This sample was incorporated into analyses of remaining replicate 393 buckets, which were specifically used to track shifts in phenotypic and genetic dynamics 394 throughout the remainder of the larval period in each treatment. 40 in all buckets). Treatment water was removed, and culture buckets were washed three times 402 with FSW to remove unsettled larvae. Individuals that remained attached to the walls of the 403 bucket were frozen and stored at -80 o C for DNA analysis. 404

Culture system and seawater chemistry 405
Larvae were reared in a temperature-controlled sea table (17.2 o C) and 0.35 µm filtered 406 and UV-sterilized seawater (FSW), pumped from 5 m depth in the bay of Villefranche. Two 407 culture systems were used consecutively to rear the larvae, both of which utilized the additions of 408 pure CO 2 gas for acidification of FSW. First, from days 0-26 the larvae were kept in a flow-409 through seawater pH-manipulation system described in Kapsenberg et al. (2017) 55 . Briefly, 410 seawater pH (pH T 8.05 and pH T 7.4) was controlled in four header tanks using a glass pH 411 electrode feedback system (IKS aquastar) and pure CO 2 gas addition and constant CO 2 -free air 412 aeration. Two header tanks were used per treatment to account for potential header tank effects. 413 Each header tank supplied water to three replicate culture buckets (drip rate of 2 L h -1 ), fitted 414 with a motorized paddle and Honeywell Durafet pH sensors for treatment monitoring (see 415

Kapsenberg et al 2017 for calibration methods). 416
On day 27 of the experiment, the flow-through system was stopped due to logistical 417 constraints and treatment conditions were maintained, in the same culture buckets, using water 418 changes every other day. For water changes 5 L of treatment seawater (70% of total volume) was 419 replaced in each culture using FSW pre-adjusted to the desired pH treatment. Seawater pH in 420 each culture bucket was measured daily, and before and after each water change. Titrando 888 56 . Accuracy of A T measurements was determined using comparison to a certified 427 reference material (Batch #151, A. Dickson, Scripps Institution of Oceanography) and ranged 428 between -0.87 and 5.3 µmol kg -1 , while precision was 1.23 µmol kg -1 (based on replicated 429 samples, n = 21). Aragonite saturation and pCO2 were calculated using pH and A T measurements 430 and the seacarb package in R 57 with dissociation constants K 1 and K 2 58 , Kf 59 and Ks 60 . Seawater 431 chemistry results are presented in the electronic supplementary, Tables S1-S2. 432

Shell Size Analysis 433
Shell size was determined as the maximum shell length parallel to the hinge using 434 brightfield microscopy and image analysis in ImageJ software. All statistical analyses were 435 conducted in R (v. 3.5.3) 61 . As larval shell length data did not pass normality tests (Shapiro-Wilk 436 test), shell-size was log-transformed to allow parametric statistical analysis. We tested the effect 437 of day, treatment, and the interaction of the two using linear mixed effects models, with day and 438 treatment as fixed effects and replicate bucket as a random effect (lmer). Effects of treatment and 439 size class on log-transformed shell length from size-separated larvae were also analyzed using a 440 linear-mixed effect model. Size class, treatment, and their interaction were fixed effects, while 441 larval bucket was a random effect. 442

DNA extraction and exome sequencing 443
We implemented exome capture, a reduced-representation sequencing approach, to MEGABLASTed to the Mytilus galloprovincialis draft genome contigs available at NCBI 456 identification tests, and a statistical test of genomic differentiation (F ST ). We visualized patterns 503 of genetic variation throughout the experiment with PCA (prcomp function in R). Prior to PCA, 504 the allele frequency matrix (the generation of which is provided in the previous section) was 505 centered and scaled using the scale function in R. Only larval samples that encompassed the full 506 phenotypic distribution within a particular bucket were included in this analysis. In other words, 507 the rows of the allele frequency matrix corresponding to larval samples that were selectively 508 segregated based on shell size were removed, and PCA was run using the day 0 larval population 509 and larval samples collected from each treatment on days 6, 26, and 43. A separate PCA was 510 then implemented using data from the day 0 larval population and all day 6 larval samples, 511 which included discrete size groups from each treatment. This analysis thus explicitly examined 512 a genomic signature of the individuals that were phenotypically distinct. 513 We next sought to identify the presence, number, and treatment-level overlap of genetic 514 variants that significantly changed in frequency between larval samples. Specifically, Fisher's 515 Exact test (FET) and the Cochran-Mantel-Haenszel (CMH) test were used to generate 516 probabilities of observed allele frequency changes, using the package Popoolation 67 in R. P-517 values for each SNP were converted to q-values in the R package qvalue 68 , and outliers were 518 identified as those SNPs with a q-value <0.01. The FET was used to identify outliers between the 519 day 0 larval population and the day 6 larval populations in each treatment (no treatment 520 replicates were available for this comparison). For all remaining allele frequency comparisons 521 (in which replicate bucket information was available), the CMH test was used to identify 522 consistent allele frequency shifts among replicate buckets. We used this test to identify 523 43. Equal densities within each treatment were assumed as no discernable difference in mortality 543 between treatments was observed (MCB, pers. obs.). 544

Gene Identification/Ontologies 545
We next sought to explore the biological pathways that were associated with survivorship 546 in each pH treatment and/or size group during the experiment. To accomplish this, we indexed 547 the outlier loci that contained variants with significant frequency changes using the annotated 548 transcriptome provided in Moreira et al. (2015). Their annotation utilized NCBI's nucleotide and 549 non-redundant, Swissprot, KEGG, and COG databases, thus providing a thorough survey of 550 potential genes and pathways associated with our candidate SNPs. We generated gene lists for 551 pH-specific loci, which were identified as loci that showed signatures of selection on all 552 sampling days and were unique to each environment. We also generated a candidate gene list for 553 loci that exhibited shared signatures of selection in each treatment. These lists thus only contain 554 robust candidate loci (loci identified as outliers in multiple independent replicates), with 555 potentially strong effect sizes (loci identified as outliers at multiple developmental stages). 556 Lastly, we used the Moreira et al. (2015) annotation to explore the genes that exhibited 557 signatures of selection for shell growth in ambient and low pH conditions, as well as shared 558 signatures of selection for shell size in each treatment. 559 560 Acknowledgements: We thank Samir Alliouane for extensive technical assistance during the 561 completion of the experiment. We would also like to thank Angelica Miglioli for experimental 562 assistance, Régis Lasbleiz for microalgal supply for larval feeding, Jacob Enk at Arbor 563 Biosciences for guidance during exome capture, and the University of Chicago Genomics Core 564 Facility for sequencing assistance. Lastly, we thank D. Rice and T. Price for insightful comments 565 on the manuscript.    Table S1. Carbonate chemistry for flow-through experimental system used on days 0-26, during which each of four 782 header tanks (two per treatment) each distributed pH adjusted seawater to three replicate buckets. Low/Amb_1/2 783 correspond to treatment replicates drawing water from separate header tanks, thus one replicate bucket per header 784 tank is represented in the table. Time-series pH and temperature data were generated using in autonomous sensor 785 were generated in each representative replicate bucket, and averge values (+/-SD) are presented . Aragonite 786 saturation (Ω a ) and pCO 2 were computed using average pH and AT for each representative replicate. Alkalinity (AT) 787 and salinity samples were generated from discrete samples taken from each of the four header tanks every other day.