Larval cultures

Mature M. galloprovincialis individuals were collected in September 2017 from the underside of a floating dock in Thau Lagoon (43.415 °N, 3.688 °E), located in Séte, France. Thau Lagoon has a mean depth of 4 m and connection to the Mediterranean Sea by three narrow channels. pH variability at the collection site during spawning season ranges from pH T 7.80 to 8.1021. Mussels were transported to the Laboratoire d’Océanographie (LOV) in Villefranche-sur-Mer, France and stored in a flow-through seawater system maintained at 15.2 °C until spawning was induced.

Within 3 weeks of the adult mussel collection, individuals were cleaned of all epibiota using a metal brush, byssal threads were cut, and mussels were warmed in seawater heated to 27 °C (~+12 °C of holding conditions) to induce spawning. Individuals that began showing signs of spawning were immediately isolated, and allowed to spawn in discrete vessels, which were periodically rinsed to remove any potential gamete contamination. Gametes were examined for viability and stored on ice (sperm) or at 16 °C (eggs). In total, gametes from 12 females and 16 males were isolated to generate a genetically diverse starting larval population. To produce pairwise crosses, 150,000 eggs from each female were placed into 16 separate vessels, corresponding to the 16 founding males. Sperm from each male was then used to fertilize the eggs in the corresponding vessel, thus eliminating the potential effects of sperm competition and ensuring that every male fertilized each female’s eggs. After at least 90% of the eggs had progressed to a four-cell stage, equal volumes from each vessel were pooled to generate the day 0 larval population (~2 million individuals), from which the replicate culture buckets were seeded. A total of 100,000 individuals were added to each culture buckets (N = 12, 18 embryos mL−1). The remaining embryos were frozen in liquid nitrogen, and stored at −80 °C for DNA analysis of the day 0 larval population. Likewise, gill tissue was collected from all founding individuals and similarly stored for downstream DNA analyses. Larvae were reared at 17.2 °C for 43 days. Starting on day 4, larvae were fed 1.6 × 108 cells of Tisochrysis lutea daily. Beginning on day 23, to account for growth and supplement diet, larvae diet was complemented with 0.2 μL of 1800 Shellfish Diet (Reed Mariculture) (days 23–28 and day 38) and ~1.6 × 108 cells of Chaetoceros gracilis (days 29–37 and 39–41). Algae were added as a pulse to each experimental replicate, twice daily (early morning, late afternoon). As the system was flow-through on days 0–26, the density of algae declined between consecutive pulses. The number of algae in each pulse was determined daily within the algal stock solution, and values reported above are the average of daily algal additions. Evidence of food consumption by larvae was indicated by observed food in larval guts, as well as substantial growth of larvae throughout the experiment.

Larval sampling

We strategically sampled larvae throughout the experiment to observe phenotypic and genetic dynamics across key developmental events, including the trochophore to D-veliger transition (day 6), the shell growth period (days 4–26), and the metamorphosis from D-veligers to settlement (days 40–43). On day 6 of the experiment, larvae were sampled from three of the six replicate buckets per treatment. A subset of larvae (N = 91–172) from each bucket was isolated to obtain shell length distributions of larvae reared in the two treatments. The remaining larvae were separated by shell size using a series of six Nitex mesh filters (70, 65, 60, 55, 50, and 20 μm; Supplementary Fig. 2b) and frozen at −80 °C. The smallest size group contained larvae arrested at the trochophore stage, and therefore unlikely to survive. The remaining five size classes isolated D-veligers from the smallest to the largest size. The shell length distribution of the larvae was used to inform, a posteriori, which combination of size classes would produce groups of the top 20% and bottom 80% of shell growers from each treatment. The relevant size groups from two replicates per treatment were then pooled for downstream DNA analysis of each phenotypic group. For the third replicate, a posteriori, all size groups were pooled in order to compute the allele frequency distribution from the entire larval population in each treatment on day 6. This sample was incorporated into analyses of remaining replicate buckets, which were specifically used to track shifts in phenotypic and genetic dynamics throughout the remainder of the larval period in each treatment.

Following size separation on day 6, the remaining replicate buckets (N = 3 per treatment) were utilized to track changing phenotypic and allele frequency distributions in the larval population through settlement. Larvae were sampled for size measurements on day 3 (N = 30–36 individuals), day 7 (N = 38–71 individuals), day 14 (N = 37–104 individuals), and day 26 (N = 49–112 individuals). Also on day 26, an additional ~1,000 larvae per replicate were frozen and stored at −80 °C pending DNA analysis. Finally, on day 43, settled individuals were sampled from each bucket (settlement was first observed on day 40 in all buckets). Treatment water was removed, and culture buckets were washed three times with FSW to remove unsettled larvae. Individuals that remained attached to the walls of the bucket were frozen and stored at −80 °C for DNA analysis.

Culture system and seawater chemistry

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

On day 27 of the experiment, the flow-through system was stopped due to logistical constraints and treatment conditions were maintained, in the same culture buckets, using water changes every other day. For water changes 5 L of treatment seawater (70% of total volume) was replaced in each culture using FSW preadjusted to the desired pH treatment. Seawater pH in each culture bucket was measured daily, and before and after each water change.

All pH measurements (calibration of Durafets used from days 0 to 26 and monitoring of static cultures from days 27 to 43) were conducted using the spectrophotometric method and purified m-cresol dye and reported on the total scale (pH T )53. Samples for total alkalinity (A T ) and salinity were taken from the header tanks every 2–3 days from days 0 to 26 and daily during the remainder of the experiment. A T was measured using an open cell titration on Metrohhm Titrando 88853. Accuracy of A T measurements was determined using comparison to a certified reference material (Batch #151, A. Dickson, Scripps Institution of Oceanography) and ranged between −0.87 and 5.3 μmol kg−1, while precision was 1.23 μmol kg−1 (based on replicated samples, n = 21). Aragonite saturation and pCO2 were calculated using pH and A T measurements and the seacarb package54 in R with dissociation constants K 1 and K 2 55, Kf56 and Ks57. Seawater chemistry results are presented in Supplementary Tables 1 and 2.

Shell-size analysis

Shell size was determined as the maximum shell length parallel to the hinge using brightfield microscopy and image analysis in ImageJ software. All statistical analyses were conducted in R (v. 3.5.3). As larval shell length data did not pass normality tests (Shapiro–Wilk test), shell-size was log-transformed to allow parametric statistical analysis. We tested the effect of day, treatment, and the interaction of the two using linear-mixed effects models, with day and treatment as fixed effects and replicate bucket as a random effect (lmer). Effects of treatment and size class on log-transformed shell length from size-separated larvae were also analyzed using a linear-mixed effect model in which size class, treatment, and their interaction were fixed effects, while larval bucket was a random effect. Significance of the fixed effects were tested against a null model using a likelihood ratio test.

DNA extraction and exome sequencing

We implemented exome capture, a reduced-representation sequencing approach, to identify SNPs and their frequency dynamics throughout the course of the experiment. Exome capture targets the protein-coding region of the genome, and thus increases the likelihood that identified polymorphisms are in or near functional loci58. Genomic DNA from each founding individual and larval sample was extracted using the EZNA Mollusc Extraction Kit, according to manufacturer’s protocol. DNA was quantified with a Qubit, and quality was determined using agarose gel, Nanodrop (260/280), and TapeStation analysis.

Genomic DNA was hybridized to a customized exome capture array designed and manufactured by Arbor Biosciences (Ann Arbor, Michigan) and using the species transcriptome provided in Moreira et al.23. Specifically, in order to design a bait set appropriate for capture of genomic DNA fragments, 90-nucleotide probe candidates were tilled every 20 nucleotides across the target transcriptome contigs. These densely tiled candidates were MEGABLASTed to the M. galloprovincialis draft genome contigs available at NCBI (GCAA_001676915.1_ASM167691v1_genomic.fna), which winnowed the candidate list to only baits with detected hits of 80 nucleotides or longer. After predicting the hybrid melting temperatures for each near-full-length hit, baits were further winnowed to those with at most two hybrids of 60 °C or greater estimated melting temperature in the M. galloprovincialis genome. This collection of highly specific baits with near-full-length hits to the draft genome were then down-sampled to a density of roughly one bait per 1.9 kbp of the final potential target space, in order to broadly sample the target while still fitting within our desired number of myBaits kit oligo limit. The final bait set comprises 100,087 oligo sequences, targeting 94,668 of the original 121,572 transcriptome contigs.

Genomic DNA from each sample was subject to standard mass estimation quality control, followed by sonication using a QSonica QR800 instrument and SPRI-based dual size selection to a target modal fragment length of 350 nucleotides. Following quantification, 300 ng total genomic DNA was taken to library preparation using standard Illumina Truseq-style end repair and adapter ligation chemistry, followed by six cycles of indexing amplification using unique eight nucleotide dual index primer pairs. For target enrichment with the custom myBaits kit, 100 ng of each founder-derived library were combined into two pools of 14 libraries each, whereas 450 ng of each embryonic and larval-pool derived library were used in individual reactions. After drying the pools or individual samples using vacuum centrifugation to 7 μL each, Arbor followed the myBaits procedure (v. 4) using the default conditions and overnight incubation to enrich the libraries using the custom probe set. After reaction cleanup, half (15 μL) of each bead-bound enriched library was taken to standard library amplification for ten cycles using Kapa HiFi polymerase. Following reaction cleanup with SPRI, each enriched library or library pool was quantified using qPCR, indicating yields between 30 and 254 ng each.

The captured libraries were sequenced at the University of Chicago Genomics Core Facility on three lanes of Illumina HiSeq 4000 using 150-bp, paired-end reads. The captured adult libraries were sequenced on an individual lane, while the 22, pooled larval samples were split randomly between the remaining two lanes. Average coverage for founding individuals was 40×, while average coverage in pooled samples was 100×.

Read trimming and variant calling

Raw DNA reads were filtered and trimmed using Trimmomatic59, and aligned to the species reference transcriptome provided in Moreira et al.23 using bowtie260. Variants in the founding individuals were identified using the Genome Analysis Toolkit’s Unified Genotyper61. These variants were filtered using VCFTOOLs62 with the following specifications: Minor Allele Frequency of 0.05, Minimum Depth of 10x, and a Maximum Variant Missing of 0.75. The resulting.vcf files provided a list of candidate bi-allelic polymorphisms to track at each time point, treatment, and phenotypic group in the larval samples. Accordingly, GATK’s Haplotype Caller was used to identify these candidate polymorphisms within each larval alignment file, and the resulting.vcf was filtered using VCFTOOLs and the following specifications: Minor Allele Frequency of 0.01, Minimum Depth of 50x, and Maximum Depth of 450x. Only variants that passed quality filtering and were identified in all larval samples (i.e., each day, treatment, and phenotypic group) were retained for downstream analyses. This process resulted in a candidate SNP list of 29,400 variants. Allele frequencies for each variant were computed as the alternate allelic depth divided by total coverage at the locus.

Allele frequency analysis

To explore how the allele frequency of the 29,400 SNPs changed in each environment throughout the course of the experiment, we used a combination of PCA, outlier loci identification tests, and a statistical test of genomic differentiation (F ST ). We visualized patterns of genetic variation throughout the experiment with PCA (prcomp function in R). Prior to PCA, the allele frequency matrix was centered and scaled using the scale function in R. Only larval samples that encompassed the full phenotypic distribution within a particular bucket were included in this analysis. In other words, the rows of the allele frequency matrix corresponding to larval samples that were selectively segregated based on shell size were removed, and PCA was run using the day 0 larval population and larval samples collected from each treatment on days 6, 26, and 43. A separate PCA was then implemented using allele frequency data from the day 0 larval population and all day 6 larval samples, which included discrete size groups from each treatment. This analysis thus explicitly examined a genomic signature of the individuals that were phenotypically distinct.

We next sought to identify the presence, number, and treatment-level overlap of genetic variants that significantly changed in frequency between larval samples. Specifically, Fisher’s Exact Test (FET) and the CMH test were used to generate probabilities of observed allele frequency changes, using the package Popoolation63 in R. P values for each SNP were converted to q values in the R package q value64, and significant SNPs were identified as those SNPs with a q value < 0.01. As the CMH test computes probabilities for SNPs based on consistent changes among replicates, it is a powerful approach to identifying significantly changed SNPs when treatment replicates are available63. Accordingly, this test was used identify significant allele frequency changes between the day 0 larval population and the day 26 ambient and low-pH treatment replicates (N = 3), the day 0 larval population and the settled individuals from ambient (N = 2) and low-pH treatment replicates (N = 3), and between the top 20% and bottom 80% of growers in each treatment (N = 2). This test thus produced a single, significant SNP list for each treatment on days 26 and 43, as well as the size-separated day 6 larvae. As treatment replicates were not available the ambient and low-pH larval population samples on day 6, an alternate contingency table test, FET, was used to obtain a list of significant SNPs for this day of sampling (with identical rank-based approach/multiple testing corrections as those used for the CMH tests). The resulting lists of significant SNPs for days 6, 26, and 43 were then compared with identify outlier loci for each treatment. Specifically, outliers were identified for each treatment as those loci containing significant SNPs on each sampling day (i.e., SNPs overlapping among a treatment’s three significant SNP lists). As the buckets sampled on day 6 were independent cultures from those buckets sampled on days 26 and 43, this process leverages both multiple independent larval cultures and sampling days to obtain a robust outlier list for each treatment.

To provide a third, independent metric of genomic change in the larval population throughout the experiment, we computed the F ST statistic for a series of comparisons. Specifically, we implemented a methods-of-moments estimator of F ST from Pool-seq data in an analysis of variance framework, as described in Hivert et al.65 (poolFstat package in R). A global (exome-wide) F ST statistic was computed pairwise between the day 0 larval population and the day 6 ambient and low pH larvae replicate buckets, day 26 ambient and low-pH larvae replicate buckets, and settled individuals from all replicate buckets in ambient and low pH. F ST was also computed to compare differentiation between phenotypic groups (top 20% and bottom 80% of growers) on day 6. As poolFstat necessitates an approximation of population size, we parameterized the model using larval counts obtained on days 0 and 26, and estimates of larval population size for days 6 and 43. The population size in each replicate bucket on day 0 was 100,000 individuals, as larval counts were conducted on the starting larval culture and an equal volume (containing 100,000 individuals) was added to each replicate bucket. To calculate population sizes on day 26, 200 mL of seawater was extracted from each replicate bucket, from which larvae were concentrated using a 70 μm mesh filter and subsequently photographed using brightfield microscopy. This mesh size was selected based on data obtained from the size selection of larvae (Supplementary Fig. 1), and chosen to filter out D-veligers that had failed to grow, and therefore survive, beyond day 6. The number of photographed larvae on day 26 thus provided a proxy for the population size at this time point: population size was estimated to be 3685, 2090, and 1183 in the ambient replicates and 2503, 1733, and 2888 in the low-pH replicates. These estimates suggest average mortality rates of 97.7 and 97.6 % in the ambient and low-pH conditions, respectively. These estimates are in accordance with previous studies that have similarly reported substantial mortality of marine bivalve larvae reared in laboratory settings66,67,68,69. Larval counts from day 26 were used to estimate a broad range of feasible population sizes on days 6 and 43. Specifically, given the observed population sizes on day 26 (reported above), F ST was computed over large ranges of population sizes on day 6 (10,000–40,000 individuals) to encompass both linear and nonlinear declines in population size throughout the growth period. While larval counts were additionally not conducted on day 43 (all larvae sampled on this day needed to be preserved for sequencing to ensure accurate allele frequency estimates), the larvae were observable to the naked eye at this time point. This allowed for the observations that there were: (1) no treatment-specific patterns in the total number of settled larvae and (2) the number of settled larvae was greater than 100 but far less than 400 individuals at this time point. F ST was thus computed using input pool sizes between 100 and 400 individuals for day 43 samples. Results demonstrating how computed values of F ST changed according to these differences in input pool size are reported in Supplementary Table 3.

Gene identification/Ontologies

We next sought to explore the biological pathways that were associated with survivorship in each pH treatment and/or size group during the experiment. To accomplish this, we indexed our list of outlier loci using the annotated transcriptome provided in Moreira et al.23. Their annotation utilizes NCBI’s nucleotide and nonredundant, Swissprot, KEGG, and COG databases, thus providing a thorough survey of potential genes and pathways associated with our candidate SNPs. We generated gene lists for pH-specific outlier loci, which were identified as loci that showed signatures of selection on all sampling days and were unique to each environment. We also generated a candidate gene list for loci that exhibited shared signatures of selection in each treatment. These lists thus only contain robust outlier loci (i.e., containing significant SNPs in multiple independent replicates), with potentially strong effect sizes (i.e., containing significant SNPs at multiple developmental stages). Lastly, we used the Moreira et al.23 annotation to explore the genes that exhibited signatures of selection for shell growth in ambient and low-pH conditions, as well as shared signatures of selection for shell size in each treatment.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.