Molecular basis of ocean acidification sensitivity and adaptation in Mytilus galloprovincialis.

Predicting the potential for species adaption to climate change is challenged by the need to identify the physiological mechanisms that underpin species vulnerability. Here, we investigated the sensitivity to ocean acidification in marine mussels during early development, and specifically the trochophore stage. Using RNA and DNA sequencing and in situ RNA hybridization, we identified developmental processes associated with abnormal development and rapid adaptation to low pH. Trochophores exposed to low pH seawater exhibited 43 differentially expressed genes. Gene annotation and in situ hybridization of differentially expressed genes point to pH sensitivity of (1) shell field development and (2) cellular stress response. Five genes within these two processes exhibited shifts in allele frequencies indicative of a potential for rapid adaptation. This case study contributes direct evidence that protecting species' existing genetic diversity is a critical management action to facilitate species resilience to climate change.

Introduction

Anthropogenic environmental change is placing unprecedented stress on marine ecosystems. Species’ capacity to survive is determined by their ability to move to more favorable habitats, exploit phenotypic plasticity to persist in the new conditions, or evolve via selection for resistant genotypes (Hoffmann and Sgrò, 2011). While phenotypic responses to global change stressors are relatively straightforward to measure, understanding what drives such responses at a molecular level and predicting the potential for those mechanisms to evolve remain a significant challenge (Franks and Hoffmann, 2012). Specifically, such efforts generally require experimental manipulations that extend across developmental stages or generations (e.g., Lim et al., 2021). Here, using a globally distributed marine bivalve exposed to changes in seawater pH associated with ocean acidification, we attempt to identify the molecular pathways underpinning phenotypic responses and determine the capacity of those pathways to evolve.

Ocean acidification is a consequence of absorption of anthropogenic carbon dioxide (CO2) emissions by seawater and is particularly harmful to calcifying marine invertebrates, such as bivalves (Kroeker et al., 2013). Conditions of low pH, including low saturation of carbonate ions, challenges the formation of calcium carbonate shells and skeletons. Furthermore, CO2 dissolved in seawater can cross biological membranes and cause intracellular acidification (Venn et al., 2013), affecting enzyme and protein function (Casey et al., 2010). Despite a diversity of ocean acidification impacts on marine metazoans and a historic focus on calcifying species (Kroeker et al., 2013), a recent meta-analysis of transcriptomic studies revealed several common molecular responses that include, but extend far beyond, biogenic calcification (Strader et al., 2020). Still, mechanistic links between molecular and phenotypic responses, which ultimately drive population demographics (e.g., individual survival, development, and growth), remain largely unknown. Resolving these links will be critical to effectively forecast the potential of natural populations to rapidly acclimate and/or adapt as ocean acidification progress (Sunday et al., 2014).

For many marine invertebrates, and particularly for molluscs (Gazeau et al., 2013), the larval stages are the most sensitive to changes in seawater carbonate chemistry (Kroeker et al., 2013). Ocean acidification induces compounding impacts on the early development of marine mussels (Figure 1). Briefly, exposure to projected declines in seawater pH leads to delayed development and both smaller and abnormal-shaped shells in D-veliger larvae (Kapsenberg et al., 2018; Kurihara, 2008; Waldbusser et al., 2015). Reduced shell size can be attributed to reduced carbonate saturation state in the calcifying space of the leading shell edge (Melzner et al., 2017). The consequence is a strong, and largely passive, dependence of calcification on external seawater chemistry (Kapsenberg et al., 2018; Melzner et al., 2017; Waldbusser et al., 2015). In contrast, the abnormal shape of D-veligers is driven by abnormal development of the underlying soft tissue, the shell field, caused by low pH conditions during the preceding trochophore stage (Kapsenberg et al., 2018). The primary developmental process at the trochophore stage is the development of the shell field (Figure 1). The shell field is part of the larval epidermis and its cells secrete an extracellular, chitin-based, organic matrix within which calcium carbonate is precipitated to form the larval shell. The shell field develops from an invagination of ectoderm cells, followed by a pore-closure and secretion of the periostracum by a rosette of cells remaining at the surface (Kniprath, 1980, 1981). At the trochophore stage, the deeper cells in the pore then evaginate back to the surface where, through mitotic divisions, the shell field expands and secretes the organic matrix of the larval shell under the larval periostracum. Interruption or abnormal progression of the shell field evagination is likely the cause for the observed indentation in the organic matrix at the center of the shell field when larvae are reared in low pH (Kapsenberg et al., 2018). As the larval shell is calcified, this indentation persists from the trochophore to the later D-veliger stage in the form of an abnormal shell hinge (Kapsenberg et al., 2018).

Figure 1

Ocean acidification impacts on early larval development

Chronologically, ocean acidification slows development, induces abnormal shell field development (hashed period; around 30 hours post-fertilization, hpf) and mantle development (dotted period, around 40 hpf) during the trochophore stage which impacts D-veliger phenotype, and reduces carbonate saturation state in the calcifying space which impairs shell growth (30 hpf onward). Chronology represents development at 14°C; hpf is hours post fertilization; PDI is Prodissoconch I; PDII is Prodissoconch II; colored larval images show the organic matrix secreted by the shell field (blue) and the calcified shell (green; scale bar is 50 μm), scale bar on D-veliger phenotypes is 30 μm. Figure is modified from supplemental information of Kapsenberg et al. (2018).

Fortunately, marine mussels are highly fecund which helps maintain substantial genetic variation in natural populations (Romiguier et al., 2014). Such variation may explain the diverse sensitivity to ocean acidification among different larval families (Frieder et al., 2017; Kapsenberg et al., 2018) and could contribute to rapid adaptation (Bitter et al., 2019). For example, previous research revealed that certain larval families are unaffected and exhibit entirely normal development in low pH conditions (Kapsenberg et al., 2018). Also, low pH conditions have shown to cause substantial changes to genetic variation in a genetically diverse larval population throughout development (Bitter et al., 2019). Such evidence of standing variation for low pH tolerance is further observed across a range of bivalve species (e.g., Ginger et al., 2013; Thomsen et al., 2017; Fitzer et al., 2019), highlighting the importance of understanding mechanisms of adaptation in this sensitive group.

To better understand the molecular underpinnings of pH sensitivity in larval development and the potential for marine mussels to evolve and overcome ocean acidification stress, the aims of this study were to (i) identify genes associated with the abnormal development of the shell field in trochophores of Mytilus galloprovincialis exposed to low pH conditions, and (ii) explore whether those genes that are physiologically responsive to low pH conditions harbor genetic variation that could drive rapid adaptation to ocean acidification. First, we performed genome-wide, RNA sequencing to conduct differential gene expression analysis (‘RNA screen’ from hereon) in larvae exposed to control and low pH conditions during the trochophore stage (see ‘Exp. 3’ in Kapsenberg et al., 2018, for pH exposures and larval phenotypes). Differentially expressed genes (DEGs) were verified with qPCR, and the localization of larval expression of DEGs was identified using in situ RNA hybridization (ISH) to provide potential links between the molecular and phenotypic responses to ocean acidification exposure. We then compared our list of DEGs to previous studies assessing bivalve gene expression patterns in response to low pH, as well as an independent list of genes shown to be subject to selection for low pH tolerance in this species (reported by Bitter et al., 2019). This latter comparison (‘DNA screen’ from hereon) aimed to directly explore the likelihood for pH-sensitive biological processes to evolve in response to ocean acidification. Our results uncover a new set of candidate genes that may underpin the pH sensitivity of larval shell field development. Furthermore, we provide evidence that genetic variation at these loci could facilitate rapid adaptation of M. galloprovincialis to ocean acidification.

Results

RNA screen: Gene expression and localization in trochophores

Larval cultures exposed to low pH conditions (pHTotal(T) 7.4) during the trochophore stage exhibited 34% abnormal development, compared to 3% in larval cultures exposed to pHT 8.1 during the trochophore stage. RNA sequencing of trochophores exposed to low pH conditions revealed 38 DEGs compared to trochophores exposed to pHT 8.1 that, after verification of each DEG’s DNA sequence, correspond to 38 different protein sequences (Table S1). Upregulation in low pH was detected in 37 of the DEGs (logFC ranged from 0.32 to 1.25), including 15 genes with annotations. Upregulation was confirmed by qPCR on 11 of the DEGs, nine of which showed statistical significance (Table 1 and Figure S1).

Table 1

Low pH-responsive gene expression, localization, and allele frequency change

Experiment	Gene	qPCRp valuea	Trochophore localization	Change in mean allele frequency (% ± SE)
RNA screen (gene expression in trochophores)	Tyr3	0.001∗	shell field	–
	VwkA	0.019∗	shell field	–
	Cdc42	0.002∗	shell field	–
	Ets96b-like	0.009∗	shell field	–
	Runx	0.001∗	shell field	–
	Steap4	0.045∗	shell field	–
	Ef2-like	0.038∗	shell field	–
	Mstn	0.026∗	muscle progenitors	–
	Unknown protein 1	0.001∗	–	–
	Slc12	0.157	shell field	–
	Unknown protein 2	0.409	shell field	–
	Cbx8-like	–	shell field	–
	Unknown protein 3	–	shell field	–
	Unknown protein 4	–	shell field	–
	Tob1	–	–	–
	ArfGap	–	–	–
	Imp-L2	–	–	–
	PrdX	–	–	–
DNA screen (selection from day 0 to 26)	Chi	0.002∗	shell field	33.1 ± 2.5
	Kif13	0.013∗	shell field	24.0 ± 2
	Unk-Hsp70-like	0.009∗	shell field	9.7 ± 4.3
	Trip11-like	0.011∗	endoderm	15 ± 1.2
RNA and DNA screen	Arih1	0.047∗	not localized	28.5 ± 1.1

Annotated genes exhibiting upregulation of expression during trochophore stage (RNA screen) and allele frequency change from day 0 to 26 (DNA screen) in larvae exposed to low pH conditions

qPCR p values on genes identified by the RNA and DNA screens both originate from qPCR performed on the same larval samples as used for RNA-seq. Asterisks indicate gene exhibited a qPCR p value < 0.05.

ISH on 14 of the DEGs demonstrated that 12 are expressed in the shell field: three genes, Tyr3, VwkA, and Runx, are known or suspected to be involved in larval shell morphogenesis (Middelbeek et al., 2010; Miglioli et al., 2019); seven genes (including VwkA) show localization in the shell field for the first time; and three genes encode unknown proteins (Table 1, Figure 2). The remaining two DEGs that did not exhibit shell field expression were localized in the shell muscle progenitors (Mstn) or did not show conspicuous localization (Arih1).

Figure 2

In situ RNA hybridization of genes in the mid-trochophore stage

Genes are those that exhibit differential expression, as determined by RNA sequencing, when trochophores are exposed to pHT 7.4. Scale bar is 20 μm.

(A) Genes previously known to be expressed in the shell field.

(B) Genes not previously known to be expressed in the shell field.

(C) Unknown genes showing shell field expression.

(D) Genes not showing shell field expression.

(E) SEM image of a representative trochophore larvae; shell field is outlined (dashed). All imaged trochophore larvae were reared in pHT 8.0 and fixed in the mid-trochophore stage (29 h postfertilization at 16.5°C).

DNA screen: Genes sensitive to pH that are prone to selection

To test whether molecular pathways sensitive to low pH conditions could lead to adaptation to ocean acidification, we compared the list of DEGs from the RNA screen to genes that exhibit signatures of selection during larval development reported by (Bitter et al., 2019). This DNA screen provided 30 outlier loci that exhibited strong signatures of selection in low pH conditions and mapped to annotated genes. Only one of these outliers, Arih1, was identified as differentially expressed in trochophore larvae exposed to low pH, according to our initial RNA screen (Table 1, Figure 3A). However, due to the high statistical stringency (or batch effects), we suspected this low overlap may, in part, be driven by false negatives in the RNA screen (Pomaznoy et al., 2019). As such, we explored the list of outliers from Bitter et al. (2019) for genes that have been shown to be pH-responsive in other marine invertebrates. We then used qPCR, and the same samples as those used for the RNA screen, to show that Chi, Trip11-like, Kif13, and Unk-Hsp70-like were significantly upregulated in response to low pH in trochophores, and so likely false negatives in our RNA screen (Figure 3). ISH revealed that Chi, Kif13, and Unk-Hsp70-like are strongly expressed in the whole shell field (Table 1, Figure 3), although only Chi has previously been known to be involved in larval shell development (Li et al., 2017; Liu et al., 2020). Trip11-like is expressed in the endoderm. Tyr1 was not detected by the RNA screen and showed no evidence of differential or localized expression at the trochophore stage (Figure S2). Thus, Tyr1 does not contribute to trochophore development, but likely instead supports larval survival in low pH conditions at a later developmental stage (Bitter et al., 2019).

Figure 3

Allele frequency changes and differential expression in genes identified in the DNA screen

(A–E) Frequency of selected allele at candidate genes on days 0 (N = 1), 6 (N = 1), and 26 (N = 3) in larvae reared in pHT 7.4. A constant of −0.28 was added to allele frequency values in panel (D) to present frequency shifts over the same range of values for all genes. All genes exhibited significant signatures of selection throughout development in low pH conditions (see STAR Methods).

(F–J) Log2-fold-expression of these genes in trochophore larvae exposed to pHT 8.1 and pHT 7.4 as determined by qPCR (faded circles indicate values individual replicate buckets, N = 6; diamonds correspond to treatment means; asterisk indicates significant shifts in expression at a p value threshold of 0.05).

(K–O) Localization of these genes based on in situ RNA hybridization in trochophore larvae exposed to pHT 8.0 (see Table 1 for localization). Scale bar is 20 μm.

The strong signatures of selection at the genomic level reported by Bitter et al. (2019) should correspond to observable changes in the frequency of specific genetic polymorphisms (alleles) throughout development in the low pH environment. Thus, to gain insight into the dynamics of low pH selection during early development, we examined changes in the frequency of single nucleotide polymorphisms (SNPs) during low pH selection from the embryo stage (day 0), through the trochophore to D-veliger transition (day 6), and during the pelagic feeding stage (day 26, Figure 3). For all genes except Unk-Hsp70-like, the most pronounced changes in allele frequency occurred between days 0 and 6 of development, which encompasses the trochophore stage, when the course for abnormal development is set (Figure 3). Specifically, the day 0 to day 6 change in allele frequency at these genes was 41% (Chit), 27% (Kif13), 28% (Trip11-like), and 32% (Arih1). These large shifts in allele frequency were maintained throughout the pelagic phase (Figure 3, Table 1).

Discussion

Predicting the impact of global environmental change on species requires understanding both an organism’s sensitivity to changing conditions, and also how that sensitivity may be mitigated through time by natural selection on standing genetic variation and an associated increase in resilience (Hoffmann and Sgrò, 2011). Here, we targeted the precise developmental point during which low pH seawater induces abnormal development in M. galloprovincialis. By performing a differential RNA-seq analysis and reviewing published literature for target genes, a total of 43 pH-responsive genes were identified in trochophore larvae. ISH performed on 18 of these genes revealed that 15 are expressed in the shell field. Of these 15 shell field genes, 13 genes were previously unknown to be expressed in the larval shell field of bivalves, which includes three genes encoding unknown proteins. Several pH-responsive genes show signatures of selection, indicating that these processes might support adaptation to ocean acidification. This includes three genes expressed in the shell field: Chi, Kif13, and Unk-Hsp70-like; and one gene expressed in the endoderm, Trip11-like, just under the shell field. While Chi is a major effector of shell morphogenesis, Unk-hsp70-like and Arih1 may mediate cellular stress response.

A core set of genes sheds light on shell field development

We previously suggested that the abnormal D-veliger phenotypes induced by low pH exposure are a result of altered development of the shell field during the preceding trochophore stage (Kapsenberg et al., 2018). Specifically, we hypothesized that abnormal shell field development is driven by abnormal or incomplete restructuring of the ectoderm during formation of the shell field. In support of this hypothesis, we found a number of DEGs that are involved in cell proliferation and migration. These include Cdc42, transcription factors Runx and Ets96b-like, and proteins Cbx8-like and E2f-like, all of which exhibit expression in the shell field. Cdc42, a Rho subfamily GTPase is known to be a central regulator of cytoskeletal organization, which supports cell shape changes and cell migration during embryonic morphogenetic movements, similar to that involved in invagination/evagination of the shell field (Murali and Rajalingam, 2014). In early embryos, Runx function is required sometimes, if not always, for both the proliferation of progenitor cells as well as the differentiation of their descendants (Coffman, 2003, 2009), while Ets96b-like is an Ets family transcription factor known to interact with Runx during gene regulation (Arman et al., 2009). The oyster ortholog of Runx (also called Runt) is also known for regulating nacre formation in pearl oysters (Zheng et al., 2017) and is expressed in the shell field of their larvae (Song et al., 2019), but was not previously known to be pH-responsive. Cbx8 and E2f proteins are part of cellular machineries regulating cell cycle progression and have been linked to proliferation of cancerous cells (Kent and Leone, 2019; Xiao et al., 2019). Upregulation of these five genes could point to the mechanisms by which ocean acidification causes the observed abnormal development of the ectoderm into the shell field (Kapsenberg et al., 2018).

Other DEGs that are expressed in the shell field appear to be involved in the development of the secreted, chitin-based organic matrix. For example, Tyr3, Chi, and VwA (von Willebrandt factor A) have previously been implicated in larval shell morphogenesis in marine bivalves: Tyr3 is a major gene involved in chitin remodeling during the formation of the organic matrix; Chi is a well-known biomineralization gene expressed in trochophores and involved in chitin metabolism (Li et al., 2017; Liu et al., 2020). VwkA is a kinase phosphorylating vwA-containing protein, which are major molluscan shell matrix proteins (Zhao et al., 2018). Our results are also consistent with a study of the oyster Crassostrea gigas, which found that the expression of Chi and Tyr is altered when larvae are reared under low pH conditions (Liu et al., 2020). All three genes are strongly expressed in the entirety of the shell field (Figures 2 and 3). Trip11-like, a thyroid receptor-interacting protein that encodes the Golgi microtubule-associated protein 210 (GMAP-210), is essential for trafficking of secreted proteins and skeletal development, and mutations in this gene have been shown to cause the rare skeletal disorder achondrogenesis type 1A (Roboti et al., 2015). Our study suggests a role for Trip11-like in shell field formation. In conclusion, differential expression of these genes suggests that the formation of the larval organic matrix is highly pH-sensitive and contributes to the delayed development associated with low pH exposure.

In addition, two proteins involved in ionic regulation were upregulated by low pH. Steap-4 (six-transmembrane epithelial antigen of the prostate 4) is a metalloreductase involved in iron and copper homeostasis (Ohgami et al., 2006). Steap-4 could play a role in the formation of the organic matrix, as iron is one of the most important minor elements in the shells of bivalves (Zhang et al., 2003; Wang et al., 2009; Coba de la Peña et al., 2016), and copper ions are necessary for tyrosinase function (Nagai et al., 2007). Slc12 is a Na+-K+-2Cl- cotransporter regulating chloride intake (Wang et al., 2009), which may be involved in the low pH responses via upregulation of HCO3/Cl pumps during high CO2 exposure in bivalves (Zhao et al., 2018). Furthermore, three unknown DEGs exhibited shell field expression and provide new candidates for future gene function research.

The RNA screen also revealed myostatin (Mstn), which is a growth factor that inhibits muscle growth in vertebrates, and has been shown to be upregulated in response to environmental changes in adult clams (Morelos et al., 2015). To our knowledge, this is the first record of larval expression of Mstn in bivalves, and its upregulation may also contribute to the observed developmental delay under low pH.

Finally, we found evidence for a cellular stress response in trochophores exposed to low pH, through the observed upregulation of Unk-Hsp70-like and Arih1. Heat shock proteins (Hsp) serve as molecular chaperones that cope with the stress-induced denaturation of other proteins (Feder and Hofmann, 1999), and shifts in either Hsp70 expression or protein production have been observed in adult M. galloprovincialis (Bitter et al., 2021), red urchins (O’Donnell et al., 2009), oysters (Timmins-Schiffman et al., 2014; Wang et al., 2016), clams (Cummings et al., 2011), sea stars (Hernroth et al., 2011), and snails (Lardies et al., 2014). However, Unk-Hsp70-like is not a bona fide Hsp70 from Mytilus, but a completely novel gene that has never been characterized and further research is needed to ascertain its roles in larval development of bivalves. Arih1 is an Ariadne1-type E3 ubiquitin ligase similar to Parkin, which binds to the substrates of misfolded or damaged proteins (Flick and Kaiser, 2012) and suppresses cell death induced by unfolded protein stress (Imai et al., 2000). Ubiquitin ligases are involved in the low pH response in other marine invertebrates including oysters (Timmins-Schiffman et al., 2014), urchins (O’Donnell et al., 2009), corals (Kaniewska et al., 2015), and coccolithophores (Jones et al., 2013). In trochophores, Arih1 exhibits expression in the endoderm, indicating that ocean acidification induces stress in non-calcifying tissues. If, and how, cellular stress might contribute to abnormal tissue development in trochophores remains unknown. It is noteworthy that in a previous study we identified differential expression of both an HSP-70 and an E3 ubiquitin ligase gene in adult M. galloprovincialis individuals (Bitter et al., 2021). Such concordance in patterns of expression indicates some degree of conserved molecular responses to low pH conditions across life stages, though more research in this area is warranted.

A core set of genes underlies adaptive potential to low pH

Accurately predicting the impacts of climate-related stressors on natural populations hinges upon identifying the extent to which adaptation may offset anticipated phenotypic consequences. Adaptation to environmental change can proceed rapidly when populations evolve via variation currently maintained within natural populations (Barrett and Schluter, 2008; Messer et al., 2016). The presence of adaptive genomic variation within M. galloprovincialis was first reported by Bitter et al. (2019), and those results have been leveraged here to assess allele frequency dynamics at specific genes with putative links to the phenotypic traits affected by low pH conditions. Based on Bitter et al. (2019), and other ocean acidification studies, we identified five genes, Chi, Trip11-like, Kif13, Unk-Hsp70-like, and Arih1, that are pH-responsive and also exhibit substantial scope for rapid adaptation in M. galloprovincialis. These genes are involved in the development of the shell field (Chi, Trip11-like, and Kif13) and the cellular stress response (Unk-Hsp70-like and Arih1).

Changes in allele frequencies at these genes from the embryo stage to the shell growth period (day 26) were substantial, ranging between 10% and 33%, indicating a strong selective advantage of those individuals harboring the adaptive allele in low pH conditions. Interestingly, for four of these genes (Chi, Trip11-like, Kif13, and Arih1), the selected allele increased most dramatically between days 0 and 6, the period encompassing the trochophore stage. This verifies that developmental processes occurring within this time window, such as shell field development and successful transition to the D-veliger stage, are important for survival and adaptation (Kapsenberg et al., 2018; Bitter et al., 2019).

The only gene exhibiting linear changes in allele frequency throughout the sampling period was Unk-Hsp70-like. While this gene was localized to the shell field, further research is needed to determine how it may participate in the bivalve stress response and to decipher its role in low pH adaptation. The cellular stress response is induced by low pH conditions in a range of marine taxa, including other marine invertebrates with similar life-history characteristics to M. galloprovincialis (e.g., high levels of genetic diversity, relatively short generation times) (Strader et al., 2020). Thus, this pathway may similarly exhibit genetic variation that supports rapid adaptation in a range of putatively sensitive species. Interestingly, a recent study identified SNPs in HSP70, as well as tyrosinase and chitinase genes, to be associated with the rate of shell growth in the oyster Crassostrea virginica, potentially corroborating the involvement of these genes in low pH adaptation (Zeng and Guo, 2021).

Conclusion

Given the focus of ocean global change biology research on characterizing species sensitivity to climate-related changes, there is a need to expand research on understanding adaptation potential (Sunday et al., 2014). Through a series of independent screens and analyses, we identified a core set of genes associated with the pH sensitivity of shell field development and found evidence for rapid adaptation potential thereof. Such adaptation could mitigate the anticipated negative impacts of ocean acidification. This study contributes new evidence to our increasing knowledge on how species could develop resilience to changing conditions (Bitter et al., 2019; Brennan et al., 2019; Kapsenberg and Cyronak, 2019). Understanding how resilience could be bolstered will help inform management practices (Kapsenberg and Cyronak, 2019). For M. galloprovincialis, the results presented here suggest that protecting species’ existing genetic diversity (e.g., establishing Marine protected areas in areas with frequent low pH exposures; Roberts et al., 2017; Kapsenberg and Cyronak, 2019) is a critical management action that could maximize the potential for rapid adaptation.

Limitations of the study

While this study’s focus on an ecologically and economically important marine bivalve in response to a global change stressor augmented its applicability to the management of natural populations, there are several limitations regarding the extrapolation of the results presented here. First, the findings from this study are limited to rapid, and immediate, adaptation as only a single generation of larvae were exposed to low pH conditions. However, it may take decades before mussel populations are exposed to low pH conditions used here, and during that time, the genetic variation in natural populations may be influenced by other stressors, such as warming or pollution. Second, our sequencing strategy (exome capture) primarily targeted coding sequences of the M. galloprovincialis genome. However, it is possible that mutations in cis-regulatory regions may influence patterns of gene expression and protein dosage in a way that influences larval development and survival. While such genotype-specific patterns of expression are known to influence adaptation in a range of species (Wittkopp and Kalay, 2012), they should be investigated further in this system to identify the causal mechanisms leading to abnormal larval phenotypes in low pH seawater.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Mark C. Bitter (mcbitter@stanford.edu).

Materials availability

Primers used for qPCR and ISH are listed in the supplemental information of this publication.

Experimental model and subject details

Mussel larvae

M. galloprovincialis gravid adults were spawned and used to create larval cultures that were then reared in the lab and sampled. Two independent experiments (RNA screen and DNA screen) were performed and both these experiments are described in detail in previous publications (Kapsenberg et al., 2018 and Bitter et al., 2019, respectively). Adult mussels were collected from two different locations: RNA screen adults were collected in February 2017 from a weather buoy line in the Bay of Villefranche-sur-Mer in the Mediterranean Sea (43.682°N, 7.319°E, Kapsenberg et al., 2018); DNA screen adults were collected in September 2017 from under a floating dock in Thau Lagoon, in Séte, France (43.415 °N, 3.688 °E, Bitter et al., 2019). A third larval culture was used to generate samples for ISH, and gravid adult mussels for these cultures were also collected from Thau Lagoon, in early October 2018. As written in Kapsenberg et al. (2018): “Adult mussels were kept in a single temperature-controlled flow-through sea table (approx. 16°C) and fed three times per week a mixture of Instant Algae (Reed Mariculture,Iso 1800 and Shellfish Diet 1800) until spawning. Spawning was induced by cleaning the mussels of epibionts and warming seawater to 28°C. Sperm was kept on ice and eggs were kept at approximately 16°C until fertilization. Test fertilizations were performed to ensure gamete compatibility prior to batch fertilization. Fertilization was performed simultaneously across biological replicates. Successful fertilization was determined by greater than 90% presence of polar bodies, after which embryos were transferred to culture vessels at a density of 14 embryos mL−1” (80,000 larvae per culture vessel) for the RNA screen experiment, and 18 embryos per mL (100,000 larvae per culture vessel) for the DNA screen experiment (Bitter et al., 2019).

Method details

RNA screen: Gene expression in trochophores

M. galloprovincialis larvae were reared in a lab system under stable and variable pH treatments as described by Kapsenberg et al. (2018) (‘Exp. 3’ therein). Embryos from nine different parental pairs (i.e., 9 males and 9 females, each used in only one cross) were pooled and distributed across four treatments that target and isolate pH exposures at the trochophore stage: (i) stable pHT 8.1, (ii) variable pHT 7.4, wherein embryos started development at pHT 7.4 but pH was increased to pHT 8.1 during the trochophore stage and reduced back to pHT 7.4 during the early D-veliger stage, (iii) stable pHT 7.4, and (iv) variable pHT 8.1, wherein embryos started development at pHT 8.1 but pH was decreased to pHT 7.4 during the trochophore stage and increased back to pHT 8.1 during the early D-veliger stage. Each treatment had three replicate culture vessels for a total of twelve larval cultures. Larvae were reared at 14.3 ºC, near habitat conditions during the spawning season in the Bay (February; Kapsenberg et al., 2017). For statistical analyses, treatments were grouped according to the pH exposure at the time of sampling during the trochophore stage, regardless of whether or not the pH treatment was fluctuating or static. This was done in accordance to phenotypic outcomes of those treatments and resulted in two groups (N = 6 per group): (i) stable pHT 8.1 and variable pHT 7.4, which combined yielded 3% abnormal development, and (ii) stable pHT 7.4 and variable pHT 8.1, which combined yielded 34% abnormal development (see Figure 3 in Kapsenberg et al., 2018).

Larvae were sampled at the mid-trochophore stage, 33 h post-fertilization (hpf). Approximately 60,000 larvae per culture vessel were preserved in TRIzol® and frozen at −80°C until RNA extractions, following manufacturer’s protocol. Quality of total RNA was verified on a NanoDrop and in an 1.5% agarose gel, and total RNA samples were aliquoted and stored at −80°C. Samples were shipped on dry ice from France to the University of Chicago where total RNA quality was checked on a bioanalyzer (mean RIN score 9.5 ± 0.3, N = 12) prior to library preparation (illumina® TruSeq Stranded mRNA). All samples were sequenced on a single lane of Illumina HiSeq4000 (100 bp, paired-end reads, average output of 6.4 megabases per sample) and a separate aliquot of total RNA was used for qPCR verification of candidate genes.

qPCR verification

DEGs identified by transcriptomic analyses were validated using qPCR. qPCR was performed on total RNA samples from all 12 cultures (4 treatments, 3 replicates) produced in Kapsenberg et al. (2018). Quantity of total RNA samples previously extracted was determined via Nanodrop and quality was verified on a subset using a bioanalyzer (RIN score range 9.6–9.7, N = 4). Total RNA was converted to cDNA using the Applied Biosystems High-Capacity RNA-to-cDNA Kit, as per manufacturer instructions. qPCR reactions were performed on the Applied Biosystems Step One real-time qPCR machine, with a total reaction volume of 15 μL consisting of 7.5 μL SYBR Green PCR Master Mix (Applied Biosystems), 1 uL 10 μM primer (Table S2), 3 μL DNA, and 3 μL DI water per reaction. As with the transcriptome analysis, six replicate samples were run for each treatment. Controls lacking cDNA template were included in each reaction plate run to identify the presence of any non-target contamination. Amplifications were performed using the Applied Biosystems StepOnePlus real-time PCR system. Melting point assessment used to further determine the efficacy of each primer pair.

In situ RNA hybridization (ISH)

To identify the location of expression of pH-responsive genes within the larval body, primers for ISH were designed from sequences when possible (Table S3) as described by Miglioli et al. (2019). Larval samples for ISH originate from an experiment independent of the transcriptomic and genomic approaches. Briefly, larvae from 4 unique male-female pairs were reared separately in static cultures (pHT 8.0, 16.5°C) and sampled at 29 hpf at which time larvae were in the mid-trochophore stage. Larvae were concentrated and fixed in 4% paraformaldehyde for 2 h at room temperature, rinsed with filtered seawater (FSW) three times, rinsed with 50% methanol for 15 min, and stored in 100% methanol at −20°C. Samples were shipped on ice for preparation for ISH following methods described in Miglioli et al. (2019) (see details in supplemental information). Samples for in situ hybridization were prepared as previously described by Miglioli et al. (2019). Briefly, 29 h post-fertilization (hpf) larvae from control pH conditions were collected and fixed overnight at 4°C with 4% Paraformaldehyde in PBS, washed in the same buffer and stored in 100% methanol at −20°C. The samples were subsequently shipped on ice and processed following the ISH protocol for M. galloprovincialis larvae described in previous work (Miglioli et al., 2019). After re-hydration in MFSW, larvae were washed 2 × 5 min in 0.1% Tween-20 in PBS-PBST, incubated in 10 μg/mL Proteinase K in PBS for 25 min, washed 2 × 5 min in PBST, fixed with 4% PFA in PBS for 40 min and washed in PBST 3 × 5 min. Samples were then placed in Hybridization buffer (50% formamide, 5X Denhart’s, 1 mg/mL yeast Torula total RNA, 0.1% Tween-20) at room temperature for 10 min, then at 65°C for 2 h. Selective probes were then added in the same solution and hybridized for 16 h at 65°C. Larvae were then washed in 5 solutions: 2 × 30 min in solution 1 (50% formamide, 5X Saline sodium citrate-SCC (20X stock solution, pH 7.0), 0,1% Sodium Dodecyl Sulfate-SDS) at 65°C, 1 × 30 min in solution 2 (50% formamide, 2X SSC, 0,1% SDS) at 65°C, 1 × 20 min in solution 3 (2X SSC, 0.1% tween 20) at 65°C, 2 × 30 min in solution 4 (0.2X SSC, 0.1% Tween 20) at 65°C, PBST 5 min at RT, 1 × 60 min at room temperature in blocking solution (0.1 M Tris, pH 7.5; 0.15 M NaCl; 0.5% Boheringer Blocking Agent -BBR, Roche, Milan). Detection was performed after overnight incubation at 4°C with an alkaline phosphatase coupled anti-DIG antibody (Sigma Aldrich, Lyon, France) diluted 1:4000 in blocking solution. After 10 × 5 min washes with 0.1% Tween-20 in 1X Tris-Buffered Saline-TBST (50 mM Tris-HCl, pH 7.6; 150 mM NaCl), larvae were washed once with the Nitro Blue Tetrazolium/5-Bromo-4-Chloro-3-Indolyl-Phosphate/ (NBT/BCIP) buffer (0.1 M Tris-HCl pH 9.5, 50 mM MgCl2, 100 mM NaCl, 0.1% Tween, 2% NBT and BCIP) and the precipitation was induced with addition of NBT/BCIP in the same buffer. Blue colour development was stopped with 1 × 5 min wash in TBST with EDTA 50 mM and 3 × 5 min washes in TBST. Samples were stored in 80% glycerol at 4°C. Images were taken with a ZEISS Axio Imager 2 microscope equipped with a B&W cooled CCD camera (Princeton Instruments MicroMax). Samples were only used if cultures of single male-female pairs yielded >90% normal development in pHT 8.0.

DNA screen: Selection throughout larval development

Data from Bitter et al. (2019) were used to identify whether or not DEGs in trochophores also exhibit adaptive variation for low pH tolerance. Briefly, the study crossed 16 males to each of 12 females and reared the resulting, genetically diverse, pooled larval population in pHT 8.1and pHT 7.4 (N = 6 per treatment) (17.2°C; note this is warmer than the temperature used in the RNA screen experiment as a warmer temperature was needed for the DNA screen to ensure timely developmental progression through settlement) from the early embryo (4-cell stage) through settlement (day 43). Larval samples were collected for allele frequency estimation on days 0, 6, 26, and 43 and sequenced using exome capture, which targets the protein-coding region of the genome, facilitating the identification of genetic polymorphisms that are in or near functional regions of the genome. Outlier loci (i.e., genes under selection) were stringently identified as those genomic regions exhibiting significant shifts in allele frequency between the day 0 larval population and the larval population on days six (N = 1), twenty-six (N = 3), and forty-three (N = 3), resulting in 30 outlier loci that mapped to functionally annotated genes (see Bitter et al., 2019 for details). The reduced replication on day six is a result of the destructive nature of the sampling, and the simultaneous need to maintain replicates throughout the duration of the experiment. As the data generated on this day of sampling were not treated as independent from subsequent sampling days (i.e. outlier loci were only identified as those genes displaying significant shifts in allele frequency on days 3, 26, and 43; see DNA Screen), this reduction in replication did not skew inferred outliers, but rather made outlier identification more stringent. These 30 genes were compared to the list of DEGs generated from the RNA screen. An additional five genes from Bitter et al. (2019), which have previously been shown to be pH-responsive in other studies/species, were tested for differential expression and location of expression in trochophores by qPCR and ISH. These genes were Chi, Trip11-like, Kif13, Hsp70, and Tyr1 (Beszteri et al., 2018; Bitter et al., 2021; Cummings et al., 2011; Glazier et al., 2020; Hernroth et al., 2011; Hüning et al., 2012; Kaniewska et al., 2015; Lardies et al., 2014; Li et al., 2017; Liu et al., 2020; O’Donnell et al., 2010; Padilla-Gamiño et al., 2013; Timmins-Schiffman et al., 2014; Wang et al., 2016). Note that for the transcript annotated as Hsp70 in Bitter et al. (2019), we refer to it as unknown-Hsp70-like (unk-Hsp70-like) given its unlikeness with the three known Hsp70 in M. galloprovincialis genome upon close evaluation.

Quantification and statistical analysis

RNA screen

Raw reads were preprocessed using the publicly available expHTS pipeline to remove PhiX contamination, trim adapters and low quality read tails and polyA’s. Reads were mapped to a reference M. galloprovincialis transcriptome using SALMON (Moreira et al., 2015; Patro et al., 2017). The resulting counts for all samples were analyzed using custom R scripts and the edgeR library (Robinson et al., 2010). After filtering, a total of 48,000 transcripts were available for the identification of DEGs. Genes were identified as differentially expressed given an FDR corrected p value < 0.05. Annotation of DEGs was verified, and in some cases modified, from those provided in Moreira et al. (2015) by re-blasting sequences on NCBI. It was not possible to find potential orthologs of the several sequences that we decided to name “unknown protein”. These sequences are identified by their gene ID (Table S1) and can be found in the NCBI protein database.

Relative expression of target mRNAs was computed using the comparative CT method described in Schmittgen and Livak (2008) and is reported as the log2-transformed fold change in low pH (pHT 7.4), relative to the control (pHT 8.1) (Schmittgen and Livak, 2008). The significance of treatment on expression (p value < 0.05) of each gene was determined using a permutational multivariate analysis of variance (PERMANOVA) in R with 999 permutations.

DNA screen

Frequency of selected alleles were quantified in the starting embryo population (day 0), three days after the trochophore to D-veliger transition (day 6), and mid-way through the pelagic feeding stage (day 26) using pooled sequencing of the larval population (Figure 2) (see Bitter et al., 2019 for further details). As described in Bitter et al. (2019), probabilities of observed shifts in allele frequencies between the starting embryo population and low pH treatment on each subsequent sampling day were generate using Fisher’s Exact Test (day 6, when allele frequencies from only one replicate were obtained) and the Cochran-Mantel-Haenszel Test (days 26 and 43, when allele frequencies from three replicates were obtained). For each sampling day, significant alleles were identified as those with shifts in frequency corresponding to a q value < 0.01. Low pH outlier loci were subsequently identified as those genes with significant alleles on days 6, 26, and 43, thus suggesting strong and consistent signatures of selection throughout larval development (Bitter et al., 2019).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Mytilus galloprovincialis	Natural population from the Bay of Villefranche-sur-Mer in the Mediterranean Sea (43.682°N, 7.319°E)	none
Bacterial and virus strains
Escherichia coli DH5 alpha	ThermoFischer	Cat#18265017
Chemicals, peptides, and recombinant proteins
TRIzol™ Reagent	ThermoFisher	Cat# 15596026
pGEM®-T Easy Vector	Promega	Cat#A1360
Anti-Digoxigenin-AP, Fab fragments	Sigma-Aldrich	Cat#11093274910
DIG nucleotides	Sigma-Aldrich	Cat#11277073910
Nitro tetrazolium Blue chloride	Sigma-Aldrich	Cat#11585029001
5-Bromo-4-chloro-3-indolyl phosphate disodium salt	Sigma-Aldrich	Cat#B6149
Tween-20	Sigma-Aldrich	Cat#P1379
Triton-X 100	Sigma-Aldrich	Cat#X100
Formamide	Sigma-Aldrich	Cat#F9037
Dodecyl sulfate sodium salt	Sigma-Aldrich	Cat#11667289001
Denhardt	Sigma-Aldrich	Cat# D9905
yeast RNA/tRNA	ThermoFisher	Cat# AM7119
Saline-Sodium Citrate buffer	Sigma-Aldrich	Cat# SRE0068
Blocking reagent	Sigma-Aldrich	Cat# 11096176001
Proteinase K	Sigma-Aldrich	Cat# 70663
Enzyme SuperScript III	Invitrogen	Cat# 18080-093
SP6 RNA polymerase	Promega	Cat# P406A
T7 RNA polymerase	Promega	Cat# P407A
Critical commercial assays
RNA extraction kit	ThermoFisher	Cat# 4387406
High-Capacity RNA-to-cDNA™ Kit	ThermoFisher	Cat# 4387406
Deposited data
DNA-seq (exome capture)	NCBI	PRJNA578149
RNA-seq	NCBI	PRJNA776346
Software and algorithms
R v.4.1.2 (see manuscript for packages used)	R Foundation for Statistical Computing
Fiji-ImageJ	https://imagej.net/	Open source
Code	GitHub	https://github.com/MarkCBitter/MolecularBasis_OA_Adaptation