# Magnitude and predictability of pH fluctuations shape plastic responses to ocean acidification

M. Bitter, Lydia E Kapsenberg, K. Silliman, Jean-Pierre Gattuso, C.A. Pfister

► To cite this version:

M. Bitter, Lydia E Kapsenberg, K. Silliman, Jean-Pierre Gattuso, C.A. Pfister. Magnitude and predictability of pH fluctuations shape plastic responses to ocean acidification. The American Naturalist, 2021, 10.1086/712930. [https://doi.org/10.1086/712937](https://doi.org/10.1086/712937)

HAL Id: hal-03043977

https://hal.science/hal-03043977

Submitted on 30 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Magnitude and Predictability of pH Fluctuations
## Shape Plastic Responses to Ocean Acidification

Mark C. Bitter,¹* Lydia Kapsenberg,²,³ Katherine Silliman,⁴ Jean-Pierre Gattuso,²,⁵ and Catherine A. Pfister¹

1. Department of Ecology and Evolution, University of Chicago, Chicago, Illinois 60637; 2. Laboratoire d'Océanographie de Villefranche, Sorbonne Université, CNRS, Villefranche-sur-Mer, France; 3. Department of Marine Biology and Oceanography, Consejo Superior de Investigaciones Científicas Institute of Marine Sciences, Barcelona, Spain; 4. School of Fisheries, Aquaculture, and Aquatic Sciences, Auburn University, Auburn, Alabama 36849; 5. Institute for Sustainable Development and International Relations, Paris, France

Submitted February 26, 2020; Accepted November 17, 2020; Electronically published February 26, 2021

Online enhancements: supplemental PDF, files. Dryad data: https://doi.org/10.5061/dryad.tjq2bvxc.

**ABSTRACT:** Phenotypic plasticity is expected to facilitate the persistence of natural populations as global change progresses. The attributes of fluctuating environments that favor the evolution of plasticity have received extensive theoretical investigation, yet empirical validation of these findings is still in its infancy. Here, we combine high-resolution environmental data with a laboratory-based experiment to explore the influence of habitat pH fluctuation dynamics on the plasticity of gene expression in two populations of the Mediterranean mussel, *Mytilus galloprovincialis*. We linked differences in the magnitude and predictability of pH fluctuations in two habitats to population-specific gene expression profiles in ambient and stressful pH treatments. Our results demonstrate population-based differentiation in gene expression plasticity, whereby mussels native to a habitat exhibiting a large magnitude of pH fluctuations with low predictability display reduced phenotypic plasticity between experimentally imposed pH treatments. This work validates recent theoretical findings on evolution in fluctuating environments, suggesting that the predictability of fluctuating selection pressures may play a predominant role in shaping the phenotypic variation observed across natural populations.

**Keywords:** phenotypic plasticity, fluctuating selection, environmental predictability, ocean acidification, bivalves.

## Introduction

The ubiquity of environmental variability has motivated decades of theoretical and empirical work aimed at determining the mechanisms facilitating the persistence of natural populations in fluctuating environments (Lewontin and Cohen 1969; Felsenstein 1976). One such mechanism is phenotypic plasticity—generally defined as the ability of a single genotype to alter its phenotype in response to a change in environmental conditions (Scheiner 1993). An extensive body of theoretical research has explored how the strength of stabilizing selection, magnitude of variation around the mean environmental state, and extent to which changes in the selective environment are predictable influence the presence of phenotypic plasticity within a population. This work has broadly demonstrated the advantage of plasticity in heterogeneous environments (Botero et al. 2015; Tufto 2015) while highlighting that the degree of plasticity expressed by individuals, oftentimes modeled as the slope of the reaction norm (Via and Lande 1985), depends critically on the extent to which changes in the selective environment are predictable (Moran 1992; Gavrilets and Scheiner 1993; DeWitt et al. 1998; Jong 1999; Tufto 2000; Reed et al. 2010; Scheiner and Holt 2012; Botero et al. 2015; Ashander et al. 2016; Bonamour et al. 2019).

Environmental predictability is defined as the presence of accurate information regarding the future state of selection and results from correlations between the environment in which a trait develops and the environment in which selection occurs (Moran 1992; Scheiner and Holt 2012). These correlations can be direct via a strong signature of positive autocorrelation in an environmental time series (Ruokolainen et al. 2009), such as those associated with diurnal fluctuations in temperature. The correlations may alternatively be indirect and arise from the presence of a second environmental variable that correlates

* Corresponding author. Present address: Department of Biology, Stanford University, Stanford, California 94305; email: mcbitter@stanford.edu.

**ORCIDs:** Bitter, https://orcid.org/0000-0001-7607-2375; Kapsenberg, https://orcid.org/0000-0002-7361-9061; Silliman, https://orcid.org/0000-0001-5964-3965; Gattuso, https://orcid.org/0000-0002-4533-4114; Pfister, https://orcid.org/0000-0003-0892-637X.

Am. Nat. 2021. Vol. 197, pp. 486–501. © 2021 by The University of Chicago. 0003-0147/2021/19704-59809$15.00. This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0), which permits non-commercial reuse of the work with attribution. For commercial use, contact journalpermissions@press.uchicago.edu.
DOI: 10.1086/712930
strongly with the agent of selection, providing a cue for an impending shift in the selective environment (Reed et al. 2010). For example, aerial emergence in coastal marine habitats due to ebbing tides (the cue) acts as a harbinger of increasing temperature (the agent of selection) and has been shown to initiate a preemptive shift of thermal physiology in intertidal mussels (Connor and Gracey 2011). Theoretical models have robustly demonstrated that even amid large fluctuations in the environment, plasticity becomes increasingly costly as the predictability of the fluctuations declines. These costs arise from plasticity-induced mismatches between the population phenotype and the optimal phenotype, which increase the likelihood of population extinction (Reed et al. 2010; Scheiner and Holt 2012; Ashander et al. 2016). Such models have ultimately shown that plasticity is selected against in highly unpredictable environments, whereby plastic genotypes are replaced by genotypes exhibiting fixed phenotypes, a strategy referred to as bet hedging (Moran 1992; Gavrillets and Scheiner 1993; DeWitt et al. 1998; Botero et al. 2015; Tufto 2015; Scheiner et al. 2020).

Empirical work has provided general support for theoretical expectations, identifying elevated plasticity within populations inhabiting environments that are characterized as generally more variable than a reference (Kingsolver and Wiernasz 1991; Kingsolver and Huey 1998; Schaum et al. 2013; Kenkel and Matz 2017). However, explicit links between specific aspects of environmental variability, such as the predictability of fluctuations, and observed levels of plasticity are sparse in the literature (but see Bonamour et al. 2019). Targeted tests of theoretical predictions are largely confined to experimental evolution-based approaches (Manenti et al. 2015; Wieczynski et al. 2018; Leung et al. 2020), which must be interpreted with caution, as such approaches may not replicate patterns of plasticity observed in nature (Kellermann et al. 2015). Thus, comparing populations that are locally adapted to divergent patterns of fluctuating selection pressures may more effectively corroborate or refute theoretical predictions. This venture is of increasing concern as global change progresses, and plasticity, if present, may serve as the initial response of populations to the associated shift in mean conditions (Hoffmann and Sgrò 2011; Seebacher et al. 2015). Additionally, the dynamics of environmental fluctuations themselves may be altered as mean conditions change (Wigley et al. 1998), further warranting investigation into which aspects of a variable environment most predominantly shape phenotypic variation (Scheiner et al. 2020). Accordingly, the present study leveraged a natural gradient in environmental variability to test the influence of fluctuation magnitude and predictability on phenotypic plasticity in an ecologically important marine bivalve subject to a global change stressor.

The global decline in seawater pH, termed “ocean acidification,” is a consequence of increasing atmospheric carbon dioxide emissions and poses extensive threats to marine systems. Laboratory-based studies assessing its expected biological effects have mounted over the past two decades, demonstrating a range of largely deleterious phenotypic effects across taxa (Kroeker et al. 2013). While the majority of such studies have imposed static pH treatments in order to quantify these effects, oceanographic data have demonstrated the dynamic nature of pH within marine systems (Wootton et al. 2008; Hofmann et al. 2011; Wootton and Pfister 2012; Kapsenberg and Hofmann 2016). This variability is pronounced in coastal environments, where local abiotic (e.g., upwelling) and biotic (e.g., photosynthesis and respiration) processes can drive substantial differences in pH regimes across temporal scales of hours to weeks and spatial scales of meters to hundreds of kilometers (Hofmann et al. 2011; Wootton and Pfister 2012; Kapsenberg and Hofmann 2016; Kwiatkowski et al. 2016). The resulting differences in low pH exposure across habitats have been shown to drive differentiation in pH tolerance between some marine populations (Hofmann et al. 2014; Kapsenberg et al. 2017c; Vargas et al. 2017; Kapsenberg and Cyronak 2019). However, the extent to which these patterns of differentiation are influenced by the temporal dynamics of pH exposure, such as how predictably pH fluctuates through time, remains unexplored.

A coastal species that has both ecological and economic value and is sensitive to ocean acidification is the Mediterranean mussel, *Mytilus galloprovincialis*. Reductions in seawater pH negatively impact *M. galloprovincialis* across life-history stages (Gazeau et al. 2014; Kapsenberg et al. 2018), although recent evidence suggests that natural populations harbor substantial genetic variation to adapt to ocean acidification (Bitter et al. 2019). Along the southern coast of France, populations of *M. galloprovincialis* are distributed throughout habitats that are likely to vary appreciably in the dynamics of environmental variability. Specifically, shallow-water lagoons in the region have limited mixing with the open Mediterranean Sea, increasing the influence of localized atmospheric and biological processes on the water bodies and driving dramatic fluctuations in abiotic conditions (Plus et al. 2003). In contrast, populations persisting along the open coastline are exposed to water masses that are more reflective of the larger Mediterranean Sea, likely mitigating frequent and large deviations in oceanographic conditions. While preliminary population genetic analysis suggests high levels of genetic connectivity between coastal and lagoon populations (Quesada et al. 1995), evidence of local adaptation within high-gene-flow marine systems has recently accrued (Sanford and Kelly 2011). Here, we leveraged the distinct abiotic environments of the lagoon and coastal
habitats to explicitly test theoretical predictions regard-
ing how exposure to, and predictability of, stressful pH
conditions influence phenotypic plasticity. We utilized
transcriptome-wide changes in gene expression as a metric
of plasticity, as shifts in gene expression are a fundamental
cellular response of an organism to a change in its environ-
ment (López-Maury et al. 2008) and patterns of expression
elucidate both genetic adaptation and acclimatized differ-
ences among individuals (Hochachka and Somero 2002).
Additionally, assays of gene expression may identify shifts
in organismal physiology that are undetectable via the ex-
amination of macroscopic traits and have become instru-
mental in characterizing the physiological consequences
of global change (DeBiasse and Kelly 2015).

Our hypotheses are illustrated in figure 1, which depicts
patterns of gene expression for each population in benign
and stressful pH treatments as a reaction norm—a model
of phenotypic plasticity in which a genotype's (or popula-
tion's mean) phenotype is plotted as a function of an envi-
ronmental gradient (Via and Lande 1985). The slope of this
relationship provides a proxy for the degree of phenotypic
plasticity, and our results are discussed in light of this model
(Via and Lande 1985). In accordance with the theoretical

predictions described above, we hypothesized that the de-
creased predictability of pH fluctuations in the lagoon
would result in individuals exhibiting reduced gene expres-
sion plasticity between experimental pH treatments. We
hypothesized that this reduction in the reaction norm slope
is facilitated by an elevation in the reaction norm intercept
(Tufto 2015), a process mechanistically driven by a “prim-
ing” of the molecular phenotype resulting from frequent
exposure to stressful conditions (Pespeni et al. 2013; Hilker
et al. 2016).

Methods

Study Site Descriptions and Environmental Monitoring

The two Mediterranean populations of *Mytilus galloprovin-
cialis* used in this study inhabited the Bay of Villefranche-
sur-Mer (hereafter referred to as the “coastal site”; Envi-
ronment Observable Littoral buoy at 43.682°N, 7.319°E)
and the Thau Lagoon (hereafter referred to as the “lagoon”;
43.4158°N, 3.6888°E), locations separated by 300 km of
coastline. The coastal site is a south-facing inlet with steep ba-
thymetry, resulting in steady mixing with the Mediterranean

<figure><img src="image_2.png" /></figure>

**Figure 1:** Conceptual overview depicting how the magnitude and predictability of environmental fluctuations may shape patterns of phenotypic plasticity. The gene expression values of the lagoon (diamonds) and coastal (circles) populations are illustrated in benign and stressful pH treatments. Lines connecting gene expression values indicate the reaction norm, the slope of which indicates the expected amount of pH stress plasticity for each population. The lagoon population, which experiences an elevation in magnitude and reduction in predictability of pH fluctuations, is expected to exhibit an increase in the baseline expression of stress response genes in the benign pH treatment and a reduced magnitude of gene expression plasticity between treatments. Note that this schematic ignores direction and simply considers the total amount of change in gene expression.
Sea (De Carlo et al. 2013). The lagoon is one of a series of shallow-water lagoons along the southern coast of France, with a mean depth of approximately 4 m and intermittent mixing with the Mediterranean Sea via two narrow channels (Plus et al. 2003). Each mussel population persists subtidally, and patterns of differentiation are thus expected to be driven by the measured oceanographic variables described below.

Autonomous oceanographic monitoring of the coastal site was conducted between September 2016 and August 2017 using SeaFET pH sensors (Martz et al. 2010) and Sea-Bird Electronic (SBE25) profilers with 60-min sampling frequencies. Monitoring of the lagoon site was conducted between October 2016 and June 2017 using SeaFET pH and SBE37 sensors with 20-min sampling frequencies. Detailed information regarding sensor deployment and calibration can be found in the supplemental PDF (available online).

The hydrogen ion concentration of all pH observations was calculated to determine the mean and magnitude of pH variability at each site. Differences in the distribution of observed pH values at each site were compared using a Kolmogorov-Smirnov test, and differences in the variance of pH at each site were tested using Hartley's maximum F-ratio test ($F_{max}$; each performed in R ver. 3.6.1). The predictability of pH fluctuations at each site was assessed using a combination of time series autocorrelation, spectral, and environmental correlation analyses. The temporal autocorrelation of pH at each site was computed across hourly and daily lag intervals using the acf function in R. These intervals provided scales of predictability that are relevant to the plastic trait measured in our study species, gene expression (Lockwood et al. 2015), and indicated the relative influence of photosynthesis/respiration (24 h) and tidal cycles (14 and 28 days) on pH variability. Spectral analysis was further used to identify the predominant pH periodicity at each site (specrum function in R). Finally, correlation coefficients between pH and temperature were computed for each site to identify differences in the reliability of habitat-based cues (cor function in R). Temperature and pH can covary in coastal marine habitats (Cyronak et al. 2019), and although we lack direct evidence that such species use temperature as a cue for impending changes in pH, differences in the correlation of the measured parameters are likely indicative of differences in a suite of additional abiotic variables (e.g., oxygen) that may serve as cues.

## Animal Collection

Adult *M. galloprovincialis* individuals, ranging in shell length between 50 and 70 mm, were collected from each habitat in October 2016. Coastal population mussels were

collected via scuba diving from a subtidal mooring in 5 m of water located within 10 m of the oceanographic sensors. The 3-m difference between the sensor location and the coastal population has a negligible influence on the mussels' pH exposure relative to those values measured by the sensor. Specifically, the predominant driver of pH variability at this site is temperature, and during the period of maximum stratification, this amounts to only a 0.10-unit pH difference between the surface and 50 m (Kapsenberg et al. 2017a). Mussels were transported dry (in netted bags) by boat to the Laboratoire d'Océanographie in Villefranche-sur-Mer (LOV), France, where they were kept in an aerated flow-through seawater tank at 16°C. Lagoon population mussels were collected by hand from the underside of a floating dock at approximately the same depth and within 10 m of oceanographic sensors. Individuals were placed in a cooler chilled with ice and transported dry by motor vehicle to LOV, where they were stored in a flow-through seawater container at 16°C. Common-garden conditioning (described below) began 3 and 7 days after the coastal and lagoon population collections, respectively.

## Experimental Design

The experiment consisted of two phases: (1) an initial common-garden conditioning of 6 weeks followed by (2) a pH treatment exposure of 5 days, after which mussels' gill tissue was sampled for transcriptomic analysis. The common-garden conditioning aimed to remove the effects of recent environmental history on individual gene expression profiles (Hochachka and Somero 2002). The 5-day exposure to stable pH treatments acclimated each individual to the experimental benign and stressful pH treatments (Thompson et al. 2012) and allowed for the quantification of individual molecular phenotypes via transcriptome-wide patterns of gene expression. The benign pH treatment corresponded to the average condition observed at both study sites throughout the oceanographic monitoring period ($pH_T = 8.1$), while the stressful pH treatment ($pH_T = 7.7$) reflected a value known to induce physiological stress on each population (Kapsenberg et al. 2018). The stressful pH treatment value was never observed at the coastal site and is exceedingly rare in the lagoon, where observations of $pH_T < 7.71$ encompassed fewer than 0.2% of all observations and never persisted for a duration of more than 6 h (fig. 2).

For common-garden conditioning, 50 adult mussels from each population were placed in a single 43.8-L flow-through seawater tank (flow rate: 7.8 L h⁻¹). To prevent mussels from aggregating, individuals from the same population were kept in 625-cm² mesh bags that were spread evenly throughout the tank (five individuals per bag). The temperature of the common garden was 16.3°C (±0.4°C)
and was aimed to mimic ambient seawater temperature in the region (table S1, available online). Mussels were fed three times per week: 15 mL of Microfeast PZ-20 (Salt Creek) for the first 2 weeks of conditioning and 8 mL of Shellfish Diet 1800 (Reed Mariculture) for the final 4 weeks. Temperature, salinity (Mettler Toledo SevenEasy Conductivity), and pH within the common garden were monitored using discrete sampling. Open-cell titration was used to determine the total alkalinity ($A_T$) necessary to calculate in situ pH (Dickson et al. 2007). Spectrophotometric measurements of pH samples were performed on discrete samples using m-cresol indicator (Dickson et al. 2007) with an accuracy of ±0.006 units $pH_T$ (determined by Kapsenberg et al. 2017b). Aragonite saturation and $pCO_2$ were computed using the pH and $A_T$ measurement and the seacarb package in R (Gattuso et al. 2016) with dissociation constants $K_1$ and $K_2$ (Lueker et al. 2000), $K_f$ (Perez and Fraga 1987), and $K_s$ (Dickson 1990; table S1).

Following common-garden conditioning, a subset of mussels from each population was randomly selected and transferred to a pH manipulation system to assess population-specific responses to the benign and stressful pH treatments. The experimental system was a flow-through seawater pH manipulation system maintained at 15.5°C and is described in detail by Kapsenberg et al. (2017a). Additional information regarding the pH manipulation system is found in the supplemental PDF. Ten individuals from each population were distributed between two independent replicate buckets per treatment ($n = 5$ individuals per bucket, $n = 10$ individuals per population per treatment) for a total of 5 days. At the end of the exposure period, individuals within each replicate bucket were removed and 0.1–0.2 mL of gill tissue was dissected and stored in 1 mL of RNAlater stabilizing solution. Sampling of each replicate bucket was randomized, and dissections of all individuals were completed within 1.5 h. The samples were maintained at room temperature for 24 h and then transferred to an –80°C freezer pending RNA isolation.

## RNA Isolation and Sequencing

Total RNA was isolated using the PureLink RNA Mini Kit with on-column DNase treatment according to manufacturer instructions for purification from animal tissue. Bioanalyzer analysis indicated extremely high-quality extracts with RNA integrity numbers ranging between 9.8 and 9.9. Before library preparation, the quantity of all extracts was determined using a qubit. We used 500 ng of total RNA per sample for 3' mRNA sequencing with the QuantSeq FWD Kit (Lexogen; Moll et al. 2014). The 3' RNAseq, also known as TagSeq, sequences only a single fragment for each transcript, allowing for shorter sequencing reads and lower sequencing depth than traditional RNAseq, and has been shown to provide more accurate estimates of transcript abundances than standard RNAseq (Lohman et al. 2016). Sequencing library quality was assessed using a Bioanalyzer high-sensitivity chip, and the 40 libraries were multiplexed and sequenced on a single lane of Illumina HiSeq 4000 at the University of Chicago Genomics Facility, yielding an average of 6.4 million 50-bp single-end reads per sample. Raw sequencing reads were mapped to an *M. galloprovincialis* reference transcriptome (Moreira et al. 2015) using Bowtie 2 (Langmead and Salzberg 2012) and yielding an overall alignment rate of 75.6%. A custom Perl script written by Misha Matz (available at https://github.com/z0on/tag-based_RNAseq) was used to count the number of reads mapping to each putative gene in the reference for all individuals.

## Identifying Population-Specific Patterns of Gene Expression and Plasticity

Differences in gene expression were used to (1) assess patterns of population-specific responses to each pH treatment, (2) explore the extent to which observed population differentiation could be related to habitat pH, and (3) determine whether pH plasticity differed between the populations. The R package DESeq2 and associated protocols (Love et al. 2014) were used to normalize and filter transcript counts and compute log₂ fold change (LFC) in expressions for each gene using a generalized linear model and a series of defined contrasts. For each contrast, genes were filtered such that the count of transcripts mapping to it were >10 across all samples, with >25 reads mapping in at least two individuals. We quantified population differentiation by computing LFCs between populations in the benign and stressful pH conditions, using the coastal population expression values as the reference ($n_{gene} = 12,809$ and 12,733 for comparisons in the benign and stressful pH treatments, respectively). Next, we quantified differences in expression plasticity by computing LFCs for each population between treatments ($n_{gene} = 12,927$ and 12,674 for quantification of the coastal and lagoon population response, respectively). Differentially expressed genes were identified as those exhibiting a Benjamini-Hochberg false discovery rate-corrected $P < .1$ (Love et al. 2014), and patterns of expression of these candidates were visualized using heat maps produced with the pheatmap package in R (Kolde 2012). The putative function of differentially expressed genes was explored with annotation data generated in the reference transcriptome provided in Moreira et al. (2015).

To visualize and quantify the molecular phenotypic differentiation of the populations in each treatment, we
conducted principal component analysis (PCA) on the
regularized log-transformed transcript counts data using
the plotPCA function in DESeq2. An ANOVA was then
conducted to statistically test for population-based dif-
ferences along the first two principal components in
each treatment. Next, we explored the extent to which
transcriptome-wide patterns of differentiation between
populations in the benign treatment were predominantly
driven by each population's native habitat. We first de-

fined a list of pH-responsive genes, those differentially ex-
pressed between pH treatments in either population. If
habitat pH was a predominant driver of the differentia-
tion in gene expression between populations, these pH-
responsive loci were expected to exhibit elevated divergence
between populations relative to all others genes differenti-
ated between populations in the benign treatment. Accord-
ingly, a permutation ANOVA (PERMANOVA) with 999 per-
mutations was used to compare the LFC distribution of these

<figure><img src="image_3.png" /><figcaption>**Figure 2:** A, B, Time series pH data collected at the coastal (A) and lagoon (B) habitats, with dashed lines indicating mean pH values across the observation window. C, Distribution of pH values observed during the pH-monitoring period at each site, with vertical lines indicating mean pH at the coastal (solid) and lagoon (dashed) sites. Arrows indicate experimentally imposed benign (black) and stressful (gray) pH treatments. D, Temporal autocorrelation in pH across hourly lag intervals for the coastal (circles) and lagoon (diamonds) habitat. Data underlying figure 2 have been deposited in the Dryad Digital Repository (https://doi.org/10.5061/dryad.tjq2bvxc; Bitter et al. 2020).</figcaption></figure>
pH-responsive loci (n = 29) with that of all genes differ-
entially expressed between populations in the benign
treatment (n = 231; RVAideMemoir package in R).

The total magnitude of pH stress response was com-
pared between populations to generate a proxy for the re-
action norm slope of the individuals from the coastal and
lagoon habitats (fig. 1). Only genes that passed filtering in
both populations' responses to low pH were considered
(n = 11,358). Population-based differences in the distri-
bution of LFC of these genes were assessed using a two-
sided Kolmogorov-Smirnov test in R. Subsequently, a
PERMANOVA with 999 permutations was used to iden-
tify significant differences in the magnitude of LFC be-
tween populations. This analysis was conducted across all
genes for which the absolute value of LFC was computed,
as well as independently for upregulated (n = 2,934) and
downregulated (n = 2,529) genes.

Results

Site-Specific Patterns of pH Fluctuation Dynamics

The mean pH value observed at each site was compara-
ble—the coastal population experienced an average pH<sub>T</sub>
of 8.07 (n = 8,045), and the lagoon population experi-

enced an average pH<sub>T</sub> of 8.11 (n = 16,992; fig. 2A, 2B). However, we quantified substantial differences in the magnitude and predictability of the pH fluctuations observed between sites. Specifically, we found a much greater range of observed pH values in the lagoon (7.47–8.39) compared with the coastal site (7.96–8.16). Additionally, a significant difference in the distributions of observed pH values (D = 0.459, P < .001) and variance across pH values (F<sub>max</sub> = 5.27, P < .001) was detected (fig. 2C). The temporal autocorrelation of pH was elevated in the coastal site across all hourly (fig. 2D) and daily (fig. S1; figs. S1–S6 are available online) intervals. Spectral analysis identified diel variability as the predominant signature of pH variability at each site (fig. S2). Temperature profiles differed appreciably between sites (fig. S2), and the coefficient of correlation between pH and temperature was dramatically reduced in the lagoon (R<sup>2</sup> = −0.17) relative to the coastal site (R<sup>2</sup> = −0.90), although the relationship between these abiotic variables was significant at each site (coastal: t = −180.86, P < .001; lagoon: t = −22.63, P < .001; fig. S3). This observed reduction in pH temporal autocorrelation and correlation between pH and temperature ultimately indicates a less predictable environment in the lagoon compared with the coastal habitat.

<figure><img src="image_4.png" /></figure>

**Figure 3:** Results from principal component analysis assessing transcriptome-wide differentiation in gene expression profiles between individuals native to the coastal (circles) and lagoon (diamonds) habitats in benign (n<sub>genes</sub> = 12,809; A) and stressful (n<sub>genes</sub> = 12,733; B) pH treatments. Data underlying figure 3 have been deposited in the Dryad Digital Repository (https://doi.org/10.5061/dryad.tjq2bvxc; Bitter et al. 2020)
## Population Differentiation in Benign and Stressful pH Conditions

PCA indicated that the molecular phenotypes of the coastal and lagoon populations were significantly differentiated in the benign pH treatment ($F_{1,18} = 14.71, P = .001$; fig. 3A) but not in the stressful treatment ($F_{1,18} = 0.55, P = .467$; fig. 3B). A total of 260 genes were differentially expressed between populations in the benign treatment (171 upregulated, 89 downregulated), while only 95 genes were differentially expressed between populations in the stressful pH treatment (55 upregulated, 40 downregulated). Distinct patterns of expression of differentially expressed genes are further visualized as heat maps (fig. S5). Of the 260 genes differentially expressed between populations in the benign treatment, 29 were differentially expressed in response to low-pH exposure. These pH-responsive genes exhibited significantly more differentiation between popu-

lations than those 231 genes not implicated in the low-pH response ($F_{1,258} = 13.62, P = .001$; fig. 4; supplemental files 1 and 2; supplemental files 1–6 are available online).

## Population-Based Differences in pH Plasticity

The observed LFC distribution generated from each population's pH stress response was distinct ($D = 0.059, P < .001$; fig. 5A), and the total magnitude of LFC across all genes was significantly greater in the coastal population, indicating reduced plasticity in the lagoon population ($F_{1,22,714} = 25.4, P = .001$; fig. 5B). When segregated by LFC direction, the magnitude of response of the coastal population was significantly greater for down-regulated ($F_{1,5,056} = 4.12, P = .043$) but not upregulated ($F_{1,5,866} = 2.05, P = .13$) genes (fig. 5B; supplemental files 3 and 4).

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Overlaid distributions of log₂ fold change of pH-responsive genes between populations (dark gray; mean value: solid line; $n_{genes} = 29$) and all other genes differentially expressed between populations (light gray; mean value: dashed line; $n_{genes} = 231$) in the benign pH treatment ($F_{1,2,581} = 13.62, P = .001$). Data underlying figure 4 have been deposited in the Dryad Digital Repository (https://doi.org/10.5061/dryad.tjqq2bvxc; Bitter et al. 2020)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5:** A, Overlaid distributions of transcriptome-wide changes in gene expression in the coastal (dark gray; mean value: solid line) and lagoon (light gray; mean value: dotted line) populations (*D* = 0.059, *P* < .001) in response to stressful pH conditions. B, Differences in the amount of upregulated (*F*<sub>1,5,866</sub> = 2.05, *P* = .13), downregulated (*F*<sub>1,5,056</sub> = 4.12, *P* = .043), and absolute (*F*<sub>1,22,714</sub> = 25.4, *P* = .001) gene expression change in response to the stressful pH conditions between the coastal (circles) and lagoon (diamonds) populations. Data underlying figure 5 have been deposited in the Dryad Digital Repository (https://doi.org/10.5061/dryad.tjqq2bvxc; Bitter et al. 2020)</figcaption></figure>

Each population exhibited a subset of genes that were differentially expressed in response to low-pH exposure: 57 differentially expressed genes (26 upregulated, 30 downregulated) were identified in the coastal population, and 49 differentially expressed genes (16 upregulated, 33 downregulated) were identified in the lagoon population (fig. S6). It is important to note that the observed responses of each population between treatments may be driven by changes in the seawater pH (the reported variable) and/or associated changes in additional biologically relevant parameters of the carbonate chemistry system (e.g., $PCO_2$ and the aragonite saturation state; Waldbusser and Salisbury 2014; Waldbusser et al. 2015a). In light of this, we explore how existing functional annotation information may lend insight into the relative contribution of these parameters to observed biological responses in “Discussion.”

## Discussion

Environmental variability plays a dominant role in shaping the ecology and evolution of natural populations. This study conducted high frequency monitoring of seawater pH to document the fine-scale and dynamic nature of pH fluctuations within coastal ecosystems and explored the extent to which it drives divergence and

shapes unique patterns of phenotypic plasticity between populations. The results presented both corroborate emerging theory and elucidate the mechanisms underlying species resilience to future climate change.

### Local Differences in pH Variability within the Northwest Mediterranean Sea

While a consistent decline in mean pH has been documented in the northwest Mediterranean Sea (Kapsenberg et al. 2017b), our high frequency monitoring of a lagoon and open coastal location demonstrates the dynamic nature of coastal carbonate chemistry across a small portion of a species biogeographic range. The dramatic pH fluctuations observed within the lagoon are driven by the shallow depth and limited mixing with the Mediterranean Sea, which act to increase the relative influence of local weather patterns and biological processes on seawater conditions (Plus et al. 2003; Waldbusser and Salisbury 2014). These processes cause the lagoon population to more frequently encounter stressful pH conditions, the occurrences of which are less predictable than the pH fluctuations observed at the coastal site. This unpredictability arises from a reduction in the autocorrelation of pH across both hourly and daily intervals, two timescales that are relevant to changes in gene expression for *Mytilus*
mussels (Lockwood et al. 2015), as well as from a reduc-
tion in the correlation of pH with temperature. Although
we report a slightly higher mean pH for the lagoon habi-
tat, this may result from the absence of time series data
during the midsummer months, the period during which
the coastal site exhibited its lowest pH values (fig. 2A). Ul-
timately, these data demonstrate that natural populations
of the species currently persist across similar mean pH
environments that differ dramatically in the magnitude
of variation. While documentation of differences in pH
variability at similar spatial scales has been previously
reported in coastal environments (Wootton and Pfister
2012; Kapsenberg and Hofmann 2016), this is the first
study to explicitly quantify differences in the predictabil-
ity of pH fluctuations occurring on a timescale that is bi-
ologically relevant for some coastal marine species and
life stages (Kapsenberg et al. 2018; Kapsenberg and
Cyronak 2019). As climate change will alter the dynamics
of environmental variability (Wigley et al. 1998; Bindoff
et al. 2019), there is increasing urgency to understand
whether and how these dynamics shape the response of
species to global change stressors.

pH Fluctuation Dynamics Shape Patterns
of Transcriptome-Wide Gene Expression Plasticity

Analyses of transcriptome-wide patterns of expression indicated that the differences in pH fluctuation dynamics observed between the coastal and lagoon habitats potentially shape patterns of plasticity and molecular phenotypic differentiation between mussel populations. By leveraging data across all genes, we observed that expression plasticity in response to low pH exposure was significantly greater in the coastal population, indicating a relative depression of the reaction norm slope in the lagoon population (fig. 1). This finding corroborates a substantial body of theoretical work, which has robustly demonstrated the fitness costs associated with maintaining plasticity in unpredictably fluctuating environments (Moran 1992; Gavrilets and Scheiner 1993; Tufto 2000; Botero et al. 2015; Ashander et al. 2016; Bonamour et al. 2019). It is important to note that these differences in plasticity are attributable to the contribution of thousands of genes, some with large, but many with subtle, shifts in expression. The relative contribution of individual genes with diverse effect sizes on organismal physiology is outside the scope of this study but presents an interesting area of future research.

We gained mechanistic insight into the processes driv-
ing the depressed plasticity of the lagoon individuals from
our comparison of the molecular phenotype of each pop-
ulation in the benign and stressful pH treatments. Specif-
ically, PCA indicated that the transcriptomic profiles of

the populations were significantly differentiated in the
benign but not the stressful pH treatment after 6 weeks
of common-garden conditioning. The increase in pheno-
typic similarity in the stressful pH treatment is likely driven
by a conserved molecular response for coping with this
shared stressor, as suggested by existing annotation infor-
mation and discussed in detail below. It is also possible
that this pattern arises from an elevation of baseline ex-
pression of stress response genes in lagoon individuals ac-
climated to the benign treatment, a mechanism reducing
the magnitude of response necessary to acclimate to the
stressful pH environment (Barshis et al. 2013; Pespeni
et al. 2013; Palumbi et al. 2014). This hypothesis is cor-
roborated by the increased divergence of pH-responsive
genes between populations in the benign treatment.
Furthermore, it is possible that the differentiation in
transcriptome-wide patterns of expression are indicative
of a physiological priming of the lagoon population, a pro-
cess in which memory of previous stress shapes responses
to future states of stress (Hilker et al. 2016). Such priming
mechanisms have been observed widely across groups,
from bacteria to plants and animals, and while oftentimes
referred to using different terms (e.g., “hardening” in re-
sponse to thermal stress exposure; Hofmann 1999), they
largely indicate the positive effect of prior stress exposure
on future responses (Hilker et al. 2016). Thus, the differ-
entiation of the populations may be driven by acclimatiz-
ation to each distinct pH habitat. However, as the dura-
tion of our common-garden conditioning is expected to,
at least in part, remove the effects of field acclimatization
on patterns of expression (Hochachka and Somero 2002;
Whitehead et al. 2011; Pespeni et al. 2013; Moyen et al.
2020) and reductions in seawater pH have indeed been
shown to induce strong selective pressure on the species
(Bitter et al. 2019), it is possible that observed differenti-
ation is also driven by local adaptation to each distinct pH
variability regime. Indeed, there is growing evidence of
local adaptation in marine systems subject to high levels
of gene flow, and selection on variation in gene expres-
sion plasticity is increasingly recognized as a fundamen-
tal driver of divergence between populations (Whitehead
and Crawford 2006; Whitehead et al. 2011; Pespeni et al.
2013). In concert, these data suggest that the lagoon popu-
lation may be limited in plasticity as ocean acidification
drives decline in mean seawater pH, increasing the pop-
ulation's reliance on genetic adaptation for persistence
(Kelly 2019).

The trends reported here lend to the growing body of
work documenting divergence in patterns of gene expres-
sion and physiology across marine populations spanning
gradients in environmental variation (Whitehead et al.
2011; Logan et al. 2012; Barshis et al. 2013; Pespeni
et al. 2013; Palumbi et al. 2014; Dong et al. 2015; Smolina
et al. 2016; Kenkel and Matz 2017; Maynard et al. 2018). Many of these studies have similarly identified differences in gene expression plasticity and baseline expression across populations, such as in coral populations exposed to extreme thermal fluctuations (Barshis et al. 2013), oyster populations exposed to recurrent low-salinity events (Maynard et al. 2018), and seaweed (Smolina et al. 2016) and limpet populations (Dong et al. 2015) spanning latitudinal gradients in mean temperature. However, quantitative analyses of the environmental time series from which focal populations originate are often lacking, precluding mechanistic links between the various attributes of fluctuating environments and observed patterns of phenotypic differentiation. Our study is an initial attempt to do so for a marine species and, in effect, has provided support for explicit theoretical predictions regarding the influence of fluctuation magnitude and predictability on the reaction norm slope and intercept. In a broader context, our findings illustrate that contemporary spatial and temporal variability in pH throughout the Mediterranean may work to maintain variation in stress tolerance across populations, a process that may facilitate rapid adaptation as ocean acidification progresses (Barrett and Schluter 2008; Bitter et al. 2019).

*Linking Patterns of Differential Expression
to Functional Responses to Ocean Acidification*

While patterns of expression across all genes allowed us to robustly test theoretical predictions, the small subset of differentially expressed genes may lend insight into the functional response of the species to stressful pH conditions and how each population copes with the distinct environmental conditions experienced in their native habitats. Shifts in pH and associated carbonate chemistry parameters, particularly PCO₂ and aragonite saturation state, impact a variety of biological processes, including calcification, metabolism, acid-base regulation, and cellular stress mechanisms (Strader et al. 2020). These broad impacts are highlighted in the present data set. For example, in response to low pH, each population differentially expressed genes mapping to perlucin, a key player in bivalve shell biomineralization (Blank et al. 2003; Wang et al. 2008). Perlucin has previously been implicated in the response to low pH in pteropods (Moya et al. 2016) and oysters (Goncalves et al. 2017), and its differential expression likely represents the direct impacts of a reduction in the aragonite saturation state on shell growth in this species (Waldbusser et al. 2015a, 2015b). Pervasive impacts of low pH exposure on the immune system were also observed in both populations via differential expression of genes mapping to *MgClq* (Gestal et al. 2010) and numerous GTPase genes of the immune-associated protein fam-

ily (Krücken et al. 2004). Such pervasive effects of low pH on immune system functioning are in accordance with recent work in bivalves (Bibby et al. 2008; Lannig et al. 2010; Wang et al. 2016; Liu et al. 2016; Castillo et al. 2017) and have been hypothesized to result from disruptions to energy metabolism caused by the hypercapnic conditions elicited from elevated PCO₂ (Pörtner 2008; Wang et al. 2016). Population-specific responses to low pH exposure were also observed. For example, the coastal population exhibited differential expression of a fatty acid-binding protein, suggestive of shifts in energy allocation associated with metabolic depression (Stewart et al. 1994). The lagoon population's response to low pH indicated the initiation of cellular stress response pathways via differential expression of heat shock protein 70 and an E3 ubiquitin ligase (Feder and Hofmann 1999; Imai et al. 2000).

Overlap in the genes differentially expressed between treatments and populations further highlights the predominant role of habitat pH in driving the observed differentiation of the coastal and lagoon populations. Specifically, of all genes differentially expressed in response to stressful pH conditions, 28% were differentially expressed between populations in the benign pH treatment. This overlap further suggests that alterations to the baseline expression of pH-responsive genes may underlie observed differences in plasticity between populations. Those genes differentiating the populations that exhibit no functional link to low-pH physiology likely indicate the suite of additional environmental variables differing between sites. For example, while differences in temperature between sites are reported here, the limited mixing of the lagoon with the larger Mediterranean Sea has also been shown to produce distinct patterns of oxygen and nutrient cycling (Plus et al. 2003). The interaction of these variables with seawater pH can alter responses to ocean acidification (Kroeker et al. 2013; Bednaršek et al. 2016; Ramajo et al. 2016) and should be recognized as potential drivers of the population differentiation reported in this study.

*Concluding Remarks*

While observations broadly linking environmental variability to the evolution of phenotypic plasticity have mounted over decades and across systems (Kingsolver and Wiernasz 1991; Kingsolver and Huey 1998; Schaum et al. 2013; Kenkel and Matz 2017), empirical studies have lagged behind theoretical research in exploring how various attributes of environmental variability promote or constrain plasticity. To address this gap, our study leveraged a natural gradient in environmental variability to demonstrate how less predictable environments depress phenotypic plasticity within natural populations, supporting recent theoretical studies (Botero et al. 2015;
Bonamour et al. 2019). In a global change context, the predictability of fluctuating environmental conditions may be altered alongside shifts in the mean environment (Wigley et al. 1998). Thus, current levels of plasticity observed within natural populations may not reflect the plasticity available to populations as global change progresses. Additional exploration into the generality of these findings across abiotic stressors and species is warranted to both robustly test theoretical predictions and better predict how species will respond to global change.

## Acknowledgments

We thank Samir Alliouane and the divers and boat crew of the Institut de la Mer de Villefranche for oceanographic sensor maintenance and mussel collection from the coastal site. We thank Remy Valdes, Behzad Mostajir, and Sébastien Mas from the Mediterranean Platform for Marine Ecosystem Experimental Research for their assistance with sensor deployment and mussel collection at the lagoon. We thank Angelica Miglioli for extensive experimental assistance and the University of Chicago Genomics Facility for sequencing assistance. We thank Nicholas VanKuren for insightful feedback on the manuscript. Last, we thank three anonymous reviewers for numerous insightful comments and suggestions on the manuscript. M.C.B. was supported by a National Science Foundation Graduate Research Fellowship Program grant (1746045) and a Department of Education grant (P200A150101). L.K. was supported by a National Science Foundation grant (OCE-1521597), which provided research support for this work. This research was also supported in part by the University of Chicago's France and Chicago Collaborating in the Sciences program to C.A.P. and M.C.B.

## Statement of Authorship

M.C.B. conceived and designed the experiment with input from L.K., J.-P.G., and C.A.P. M.C.B. and L.K. performed the experiment. M.C.B. and K.S. completed the molecular lab work, with sequencing assistance from the University of Chicago Genomics Facility. M.C.B. completed bioinformatics, statistical, and computational analyses with assistance from K.S. M.C.B. wrote the manuscript with input from L.K., K.S., J.-P.G., and C.A.P.

## Data and Code Availability

All raw data and referenced supplemental files in this article have been deposited in the Dryad Digital Repository (https://doi.org/10.5061/dryad.tjq2bvxc; Bitter et al. 2020). All code associated with statistical analyses and figure generation for this article are publicly available at GitHub (https://github.com/MarkCBitter/pHFluctuation_Plasticity) and Zenodo (https://doi.org/10.5281/zenodo.4306829; Bitter 2020).

## Literature Cited

Ashander, J., L.-M. Chevin, and M. L. Baskett. 2016. Predicting evolutionary rescue via evolving plasticity in stochastic environments. Proceedings of the Royal Society B 283:20161690. https://doi.org/10.1098/rspb.2016.1690.

Barrett, R. D. H., and D. Schlüter. 2008. Adaptation from standing genetic variation. Trends in Ecology and Evolution 23:38–44. https://doi.org/10.1016/j.tree.2007.09.008.

Barshis, D. J., J. T. Ladner, T. A. Oliver, F. O. Seneca, N. Traylor-Knowles, and S. R. Palumbi. 2013. Genomic basis for coral resilience to climate change. Proceedings of the National Academy of Sciences of the USA 110:1387–1392. https://doi.org/10.1073/pnas.1210224110.

Bednaršek, N., C. J. Harvey, I. C. Kaplan, R. A. Feely, and J. Možina. 2016. Pteropods on the edge: cumulative effects of ocean acidification, warming, and deoxygenation. Progress in Oceanography 145:1–24. https://doi.org/10.1016/j-pocean.2016.04.002.

Bibby, R., S. Widdicombe, H. Parry, J. Spicer, and R. Pipe. 2008. Effects of ocean acidification on the immune response of the blue mussel *Mytilus edulis*. Aquatic Biology 2:67–74. https://doi.org/10.3354/ab00037.

Bindoff, N. L., W. L. Cheung, J. G. Kairo, J. Arístegui, V. A. Guinder, R. Hallberg, N. Hilmi, et al. 2019. Changing ocean, marine ecosystems, and dependent communities. Pages 447–587 in IPCC special report on the ocean and cryosphere in a changing climate. Intergovernmental Panel on Climate Change, Geneva.

Bitter, M. C. 2020. MarkCBitter/pHFluctuation_Plasticity: Bitter et al. (2020) Am Nat code. Version v1.0.0. Zenodo, https://doi.org/10.5281/zenodo.4306829.

Bitter, M. C., L. Kapsenberg, J.-P. Gattuso, and C. A. Pfister. 2019. Standing genetic variation fuels rapid adaptation to ocean acidification. Nature Communications 10:1–10. https://doi.org/10.1038/s41467-019-13767-1.

Bitter, M. C., L. Kapsenberg, K. Silliman, J.-P. Gattuso, and C. A. Pfister. 2020. Data from: Magnitude and predictability of pH fluctuations shape plastic responses to ocean acidification. American Naturalist, Dryad Digital Repository, https://doi.org/10.5061/dryad.tjq2bvxc.

Blank, S., M. Arnoldi, S. Khoshnavaz, L. Treccani, M. Kuntz, K. Mann, G. Grathwohl, and M. Fritz. 2003. The nacre protein perlucin nucleates growth of calcium carbonate crystals. Journal of Microscopy 212:280–291. https://doi.org/10.1111/j.1365-2818.2003.01263.x.

Bonamour, S., L.-M. Chevin, A. Charmantier, and C. Teplitsky. 2019. Phenotypic plasticity in response to climate change: the importance of cue variation. Philosophical Transactions of the Royal Society B 374:20180178. https://doi.org/10.1098/rstb.2018.0178.

Botero, C. A., F. J. Weissing, J. Wright, and D. R. Rubenstein. 2015. Evolutionary tipping points in the capacity to adapt to environmental change. Proceedings of the National Academy of Sciences of the USA 112:184–189. https://doi.org/10.1073/pnas.1408589111.

Castillo, N., L. M. Saavedra, C. A. Vargas, C. Gallardo-Escárate, and C. Dérée. 2017. Ocean acidification and pathogen exposure modulate the immune response of the edible mussel *Mytilus* (https://github.com/MarkCBitter/pHFluctuation_Plasticity) and Zenodo (https://doi.org/10.5281/zenodo.4306829; Bitter 2020).
*chilensis*. Fish and Shellfish Immunology 70:149–155. https://doi.org/10.1016/j.fsi.2017.08.047.

Connor, K. M., and A. Y. Gracey. 2011. Circadian cycles are the dominant transcriptional rhythm in the intertidal mussel *Mytilus californianus*. Proceedings of the National Academy of Sciences of the USA 108:16110–16115. https://doi.org/10.1073/pnas.111076108.

Cyronak, T., Y. Takeshita, T. A. Courtney, E. H. DeCarlo, B. D. Eyre, D. I. Kline, T. R. Martz, et al. 2019. Diel temperature and pH variability scale with depth across diverse coral reef habitats. Limnology and Oceanography Letters 5:193–203. https://doi.org/10.1002/lol2.10129.

DeBiasse, M. B., and M. W. Kelly. 2015. Plastic and evolved responses to global change: what can we learn from comparative transcriptomics? Journal of Heredity 107:71–81. https://doi.org/10.1093/jhered/esv073.

De Carlo, E. H., L. Mousseau, O. Passafiume, P. S. Drupp, and J.-P. Gattuso. 2013. Carbonate chemistry and air-sea CO₂ flux in a NW Mediterranean Bay over a four-year period: 2007–2011. Aquatic Geochemistry 19:399–442. https://doi.org/10.1007/s10498-013-9217-4.

de Jong, G. 1999. Unpredictable selection in a structured population leads to local genetic differentiation in evolved reaction norms. Journal of Evolutionary Biology 12:839–851. https://doi.org/10.1046/j.1420-9101.1999.00118.x.

DeWitt, T. J., A. Sih, and D. S. Wilson. 1998. Costs and limits of phenotypic plasticity. Trends in Ecology and Evolution 13:77–81. https://doi.org/10.1016/S0169-5347(97)01274-3.

Dickson, A. G. 1990. Standard potential of the reaction: AgCl(s) + 12H₂(g) = Ag(s) + HCl(Aq), and the standard acidity constant of the ion HSO₄⁻ in synthetic sea water from 273.15 to 318.15 K. Journal of Chemical Thermodynamics 22:113–127. https://doi.org/10.1016/0021-9614(90)90074-Z.

Dickson, A. G., C. L. Sabine, and J. R. Christian. 2007. Guide to best practices for ocean CO₂ measurements. PICES Special Publication 3. North Pacific Marine Science Organization, Sidney, British Columbia.

Dong, Y., G. Han, M. Ganmanee, and J. Wang. 2015. Latitudinal variability of physiological responses to heat stress of the intertidal limpet *Cellana toreuma* along the Asian coast. Marine Ecology Progress Series 529:107–119. https://doi.org/10.3354/meps11303.

Feder, M. E., and G. E. Hofmann. 1999. Heat-shock proteins, molecular chaperones, and the stress response: evolutionary and ecological physiology. Annual Review of Physiology 61:243–282. https://doi.org/10.1146/annurev.physiol.61.1.243.

Felsenstein, J. 1976. The theoretical population genetics of variable selection and migration. Annual Review of Genetics 10:253–280. https://doi.org/10.1146/annurev.ge.10.120176.001345.

Gattuso, J.-P., J.-M. Epitalon, and H. Lavigne. 2016. Seacarb: seawater carbonate chemistry. R package version 3.2.2. https://cran.r-project.org/package/seacarb.

Gavrilets, S., and S. M. Scheiner. 1993. The genetics of phenotypic plasticity. V. Evolution of reaction norm shape. Journal of Evolutionary Biology 6:31–48. https://doi.org/10.1046/j.1420-9101.1993.6010031.x.

Gazeau, F., S. Alliouane, C. Bock, L. Bramanti, M. López Correa, M. Gentile, T. Hirse, H.-O. Pörtner, and P. Ziveri. 2014. Impact of ocean acidification and warming on the Mediterranean mussel (*Mytilus galloprovincialis*). Frontiers in Marine Science 1:62. https://doi.org/10.3389/fmars.2014.00062.

Gestal, C., A. Pallavicini, P. Venier, B. Novoa, and A. Figueras. 2010. MgClq, a novel C1q-domain-containing protein involved in the immune response of *Mytilus galloprovincialis*. Developmental and Comparative Immunology 34:926–934. https://doi.org/10.1016/j.dci.2010.02.012.

Goncalves, P., E. L. Thompson, and D. A. Raftos. 2017. Contrasting impacts of ocean acidification and warming on the molecular responses of CO₂-resilient oysters. BMC Genomics 18:431. https://doi.org/10.1186/s12864-017-3818-z.

Hilker, M., J. Schwachtje, M. Baier, S. Balazadeh, I. Bäurle, S. Geiselhardt, D. K. Hincha, et al. 2016. Priming and memory of stress responses in organisms lacking a nervous system. Biological Reviews 91:1118–1133. https://doi.org/10.1111/brv.12215.

Hochachka, P. W., and G. N. Somero. 2002. Biochemical adaptation: mechanism and process in physiological evolution. Oxford University Press, New York.

Hoffmann, A. A., and Sgrò C. M. 2011. Climate change and evolutionary adaptation. Nature 470:479–485. https://doi.org/10.1038/nature09670.

Hofmann, G. E. 1999. Ecologically relevant variation in induction and function of heat shock proteins in marine organisms. American Zoologist 39:889–900. https://doi.org/10.1093/icb/39.6.889.

Hofmann, G. E., T. G. Evans, M. W. Kelly, J. L. Padilla-Gamiño, C. A. Blanchette, L. Washburn, F. Chan, et al. 2014. Exploring local adaptation and the ocean acidification seascape—studies in the California Current Large Marine Ecosystem. Biogeosciences 11:1053–1064. https://doi.org/10.5194/bg-11-1053-2014.

Hofmann, G. E., J. E. Smith, K. S. Johnson, U. Send, L. A. Levin, F. Micheli, A. Paytan, et al. 2011. High-frequency dynamics of ocean pH: a multi-ecosystem comparison. PLoS ONE 6:e28983. https://doi.org/10.1371/journal.pone.0028983.

Imai, Y., M. Soda, and R. Takahashi. 2000. Parkin suppresses unfolded protein stress-induced cell death through its E3 ubiquitin-protein ligase activity. Journal of Biological Chemistry 275:35661–35664. https://doi.org/10.1074/jbc.C000447200.

Kapsenberg, L., S. Alliouane, F. Gazeau, L. Mousseau, and J.-P. Gattuso. 2017a. Coastal ocean acidification and increasing total alkalinity in the northwestern Mediterranean Sea. Ocean Science 13:411–426. https://doi.org/10.5194/os-13-411-2017.

Kapsenberg, L., E. E. Bockmon, P. J. Bresnahan, K. J. Kroeker, J.-P. Gattuso, and T. R. Martz. 2017b. Advancing ocean acidification biology using Durafet® pH electrodes. Frontiers in Marine Science 4:321. https://doi.org/10.3389/fmars.2017.000321.

Kapsenberg, L., and T. Cyronak. 2019. Ocean acidification refugia in variable environments. Global Change Biology 25:3201–3214. https://doi.org/10.1111/gcb.14730.

Kapsenberg, L., and G. E. Hofmann. 2016. Ocean pH time-series and drivers of variability along the northern Channel Islands, California, USA. Limnology and Oceanography 61:953–968. https://doi.org/10.1002/lno.10264.

Kapsenberg, L., A. Miglioli, M. C. Bitter, E. Tambutté, R. Dumollard, and J.-P. Gattuso. 2018. Ocean pH fluctuations affect mussel larvae at key developmental transitions. Proceedings of the Royal Society B 285:20182381. https://doi.org/10.1098/rspb.2018.2381.

Kapsenberg, L., D. K. Okamoto, J. Dutton, and G. E. Hofmann. 2017c. Sensitivity of sea urchin fertilization to pH varies across a natural pH mosaic. Ecology and Evolution 7:1737–1750. https://doi.org/10.1002/ece3.2776.

Kellermann, V., A. A. Hoffmann, T. N. Kristensen, N. N. Moghadam, and V. Loeschcke. 2015. Experimental evolution
under fluctuating thermal conditions does not reproduce patterns of adaptive clinal differentiation in *Drosophila melanogaster*. American Naturalist 186:582–593. https://doi.org/10.1086/683252.

Kelly, M. 2019. Adaptation to climate change through genetic ac-
commodation and assimilation of plastic phenotypes. Philosoph-
ical Transactions of the Royal Society B 374:20180176. https://
doi.org/10.1098/rstb.2018.0176.

Kenkel, C. D., and M. V. Matz. 2017. Gene expression plasticity as a mechanism of coral adaptation to a variable environment. Nature Ecology and Evolution 1:0014. https://doi.org/10.1038/s41559-016-0014.

Kingsolver, J. G., and R. B. Huey. 1998. Evolutionary analyses of morphological and physiological plasticity in thermally variable environments. Integrative and Comparative Biology 38:545–560. https://doi.org/10.1093/icb/38.3.545.

Kingsolver, J. G., and D. C. Wiernasz. 1991. Seasonal polyphenism in wing-melanin pattern and thermoregulatory adaptation in *Pieris* butterflies. American Naturalist 137:816–830. https://doi.org/10.1086/285195.

Kolde, R. 2012. pheatmap: pretty heatmaps. R Package. R Founda-
tion for Statistical Computing, Vienna.

Kroeker, K. J., R. L. Kordas, R. Crim, I. E. Hendriks, L. Ramajo, G. S. Singh, C. M. Duarte, and J.-P. Gattuso. 2013. Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming. Global Change Biology 19:1884–1896. https://doi.org/10.1111/gcb.12179.

Krücken, J., R. M. U. Schroetel, I. U. Müller, N. Saïdani, P. Marinovski, W. P. M. Benten, O. Stamm, and F. Wunderlich. 2004. Comparative analysis of the human *gimap* gene clustere encoding a novel GTPase family. Gene 341:291–304. https://doi.org/10.1016/j.gene.2004.07.005.

Kwiatkowski, L., B. Gaylord, T. Hill, J. Hosfelt, K. J. Kroeker, Y. Nebuchina, A. Ninokawa, et al. 2016. Nighttime dissolution in a temperate coastal ocean ecosystem increases under acidification. Scientific Reports 6:22984. https://doi.org/10.1038/srep22984.

Langmead, B., and S. L. Salzberg. 2012. Fast gapped-read align-
ment with Bowtie 2. Nature Methods 9:357–359. https://doi.org
/10.1038/nmeth.1923.

Lannig, G., S. Eilers, H. O. Pörtner, I. M. Sokolova, and C. Bock. 2010. Impact of ocean acidification on energy metabolism of oyster, *Crassostrea gigas*—changes in metabolic pathways and thermal response. Marine Drugs 8:2318–2339. https://doi.org/10.3390/md8082318.

Leung, C., M. Rescan, D. Grulois, and L.-M. Chevin. 2020. Reduced phenotypic plasticity evolves in less predictable environments. Ecology Letters 23:1664–1672. https://doi.org/10.1111/ele.13598.

Lewontin, R. C., and D. Cohen. 1969. On population growth in a randomly varying environment. Proceedings of the National Academy of Sciences of the USA 62:1056–1060.

Liu, S., W. Shi, C. Guo, X. Zhao, Y. Han, C. Peng, X. Chai, and G. Liu. 2016. Ocean acidification weakens the immune response of blood clam through hampering the NF-Kappa β and toll-like receptor pathways. Fish and Shellfish Immunology 54:322–327. https://doi.org/10.1016/j.fsi.2016.04.030.

Lockwood, B. L., K. M. Connor, and A. Y. Gracey. 2015. The en-
vironmentally tuned transcriptomes of *Mytilus* mussels. Journal
of Experimental Biology 218:1822–1833. https://doi.org/10.1242
/jeb.118190.

Logan, C. A., L. E. Kost, and G. N. Somero. 2012. Latitudinal dif-
ferences in *Mytilus californianus* thermal physiology. Marine

Ecology Progress Series 450:93-105. https://doi.org/10.3354/meps09491.

Lohman, B. K., J. N. Weber, and D. I. Bolnick. 2016. Evaluation of TagSeq, a reliable low-cost alternative for RNAseq. Molecular Ecology Resources 16:1315-1321. https://doi.org/10.1111/1755-0998.12529.

López-Maury, L., S. Marguerat, and J. Bähler. 2008. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nature Reviews Genetics 9:583–593. https://doi.org/10.1038/nrg2398.

Love, M. I., W. Huber, and S. Anders. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550. https://doi.org/10.1186/s13059-014-0550-8.

Lueker, T. J., A. G. Dickson, and C. D. Keeling. 2000. Ocean pCO₂ calculated from dissolved inorganic carbon, alkalinity, and equations for K₁ and K₂; validation based on laboratory measurements of CO₂ in gas and seawater at equilibrium. Marine Chemistry 70:105–119. https://doi.org/10.1016/S0304-4203(00)00022-0.

Manenti, T., V. Loeschcke, N. N. Moghadam, and J. G. Sørensen. 2015. Phenotypic plasticity is not affected by experimental evolution in constant, predictable or unpredictable fluctuating thermal environments. Journal of Evolutionary Biology 28:2078–2087. https://doi.org/10.1111/jeb.12735.

Martz, T. R., J. G. Connery, and K. S. Johnson. 2010. Testing the Honeywell DuraFet® for seawater pH applications. Limnology and Oceanography Methods 8:172–184. https://doi.org/10.4319/lom.2010.8.172.

Maynard, A., J. M. Bible, M. H. Pespeni, E. Sanford, and T. G. Evans. 2018. Transcriptomic responses to extreme low salinity among locally adapted populations of Olympia oyster (*Ostrea lurida*). Molecular Ecology 27:4225–4240. https://doi.org/10.1111/mec.14863.

Moll, P., M. Ante, A. Seitz, and T. Reda. 2014. QuantSeq 3' mRNA sequencing for RNA quantification. Nature Methods 11:i–iii. https://doi.org/10.1038/nmeth.f.376.

Moran, N. A. 1992. The evolutionary maintenance of alternative phenotypes. American Naturalist 139:971-989. https://doi.org/10.1086/285369.

Moreira, R., P. Pereiro, C. Canchaya, D. Posada, A. Figueras, and B. Novoa. 2015. RNA-Seq in *Mytilus galloprovincialis*: comparative transcriptomics and expression profiles among different tissues. BMC Genomics 16:728. https://doi.org/10.1186/s12864-015-1817-5.

Moya, A., E. L. Howes, T. Lacoue-Labarthe, S. Forêt, B. Hanna, M. Medina, P. L. Munday, et al. 2016. Near-future pH conditions severely impact calcification, metabolism and the nervous system in the pteropod *Heliconoides inflatus*. Global Change Biology 22:3888–3900. https://doi.org/10.1111/gcb.13350.

Moyen, N. E., G. N. Somero, and M. W. Denny. 2020. Mussel ac-
climatization to high, variable temperatures is lost slowly upon
transfer to benign conditions. Journal of Experimental Biology
223:jeb222893. https://doi.org/10.1242/jeb.222893.

Palumbi, S. R., D. J. Barshis, N. Traylor-Knowles, and R. A. Bay. 2014. Mechanisms of reef coral resistance to future climate change. Science 344:895–898. https://doi.org/10.1126/science.1251336.

Perez, F., and F. Fraga. 1987. The pH measurements in seawater on the NBS Scale. Marine Chemistry 21:315–327. https://doi.org/10.1016/0304-4203(87)90054-5.

Pespeni, M. H., B. T. Barney, and S. R. Palumbi. 2013. Differences in the regulation of growth and biomineralization genes revealed through long-term common-garden acclimation and experimental
genomics in the purple sea urchin. Evolution 67:1901-1914. https://doi.org/10.1111/evo.12036.

Plus, M., A. Chapelle, P. Lazure, I. Auby, G. Levavasseur, M. Verlaque, T. Belsher, J.-M. Deslous-Paoli, J.-M. Zaldivar, and C. N. Murray. 2003. Modeling of oxygen and nitrogen cycling as a function of macrophyte community in the Thau lagoon. Continental Shelf Research 23:1877-1898. https://doi.org/10.1016/j.csr.2003.03.001.

Pörtner, H.-O. 2008. Ecosystem effects of ocean acidification in times of ocean warming: a physiologist's view. Marine Ecology Progress Series 373:203-217. https://doi.org/10.3354/mep s07768.

Quesada, H., C. Zapata, and G. Alvarez. 1995. A multilocus allozyme discontinuity in the mussel *Mytilus galloprovincialis*: the interaction of ecological and life-history factors. Oceanographic Literature Review 9:769-770.

Ramajo, L., E. Pérez-León, I. E. Hendriks, N. Marbà, D. Krause-Jensen, M. K. Sejr, M. E. Blicher, N. A. Lagos, Y. S. Olsen, and C. M. Duarte. 2016. Food supply confers calcifiers resistance to ocean acidification. Scientific Reports 6:19374. https://doi.org/10.1038/srep19374.

Reed, T. E., R. S. Waples, D. E. Schindler, J. J. Hard, and M. T. Kinnison. 2010. Phenotypic plasticity and population viability: the importance of environmental predictability. Proceedings of the Royal Society B 277:3391-3400. https://doi.org/10.1098/rspb.2010.0771.

Ruokolainen, L., A. Lindén, V. Kaitala, and M. S. Fowler. 2009. Ecological and evolutionary dynamics under coloured environmental variation. Trends in Ecology and Evolution 24:555-563. https://doi.org/10.1016/j.tree.2009.04.009.

Sanford, E., and M. W. Kelly. 2011. Local adaptation in marine invertebrates. Annual Review of Marine Science 3:509-535. https://doi.org/10.1146/annurev-marine-120709-142756.

Schaum, E., B. Rost, A. J. Millar, and S. Collins. 2013. Variation in plastic responses of a globally distributed picoplankton species to ocean acidification. Nature Climate Change 3:298-302. https://doi.org/10.1038/nclimate1774.

Scheiner, S. M. 1993. Genetics and evolution of phenotypic plasticity. Annual Review of Ecology and Systematics 24:35-68. https://doi.org/10.1146/annurev.es.24.110193.000343.

Scheiner, S. M., M. Barfield, and R. D. Holt. 2020. The genetics of phenotypic plasticity. XVII. Response to climate change. Evolutionary Applications 13:388-399. https://doi.org/10.1111/eva.12876.

Scheiner, S. M., and R. D. Holt. 2012. The genetics of phenotypic plasticity. X. Variation versus uncertainty. Ecology and Evolution 2:751-767. https://doi.org/10.1002/ece3.217.

Seebacher, F., C. R. White, and C. E. Franklin. 2015. Physiological plasticity increases resilience of ectothermic animals to climate change. Nature Climate Change 5:61-66. https://doi.org/10.1038/nclimate2457.

Smolina, I., S. Kollias, A. Jueterbock, J. A. Coyer, and G. Hoarau. 2016. Variation in thermal stress response in two populations of the brown seaweed, *Fucus distichus*, from the Arctic and subarctic intertidal. Royal Society Open Science 3:150429. https://doi.org/10.1098/rsos.150429.

Stewart, J. M., R. C. Carlin, J. A. MacDonald, and S. Van Iderstine. 1994. Fatty acid binding proteins and fatty acid catabolism in marine invertebrates: peroxisomal β-oxidation. Invertebrate Reproduction and Development 25:73-82. https://doi.org/10.1080/07924259.1994.9672370.

Strader, M. E., J. M. Wong, and G. E. Hofmann. 2020. Ocean acidification promotes broad transcriptomic responses in marine metazoans: a literature survey. Frontiers in Zoology 17:7. https://doi.org/10.1186/s12983-020-0350-9.

Thompson, E. L., D. A. Taylor, S. V. Nair, G. Birch, R. Coleman, and D. A. Raftos. 2012. Optimal acclimation periods for oysters in laboratory-based experiments. Journal of Molluscan Studies 78:304-307. https://doi.org/10.1093/mollus/eys012.

Tufto, J. 2000. The evolution of plasticity and nonplastic spatial and temporal adaptations in the presence of imperfect environmental cues. American Naturalist 156:121-130. https://doi.org/10.1086/303381.

______. 2015. Genetic evolution, plasticity, and bet-hedging as adaptive responses to temporally autocorrelated fluctuating selection: a quantitative genetic model. Evolution 69:2034-2049. https://doi.org/10.1111/evo.12716.

Vargas, C. A., N. A. Lagos, M. A. Lardies, C. Duarte, P. H. Manríquez, V. M. Aguilera, B. Broitman, S. Widdicombe, and S. Dupont. 2017. Species-specific responses to ocean acidification should account for local adaptation and adaptive plasticity. Nature Ecology and Evolution 1:0084. https://doi.org/10.1038/s41559-017-0084.

Via, S., and R. Lande. 1985. Genotype-environment interaction and the evolution of phenotypic plasticity. Evolution 39:505-522. https://doi.org/10.1111/j.1558-5646.1985.tb00391.x.

Waldbusser, G. G., B. Hales, C. J. Langdon, B. A. Haley, P. Shrader, E. L. Bruner, M. W. Gray, C. A. Miller, and I. Gimenez. 2015a. Saturation-state sensitivity of marine bivalve larvae to ocean acidification. Nature Climate Change 5:273-280. https://doi.org/10.1038/nclimate2479.

Waldbusser, G. G., B. Hales, C. J. Langdon, B. A. Haley, P. Schrader, E. L. Brunner, M. W. Gray, C. A. Miller, I. Gimenez, and G. Hutchinson. 2015b. Ocean acidification has multiple modes of action on bivalve larvae. PLoS ONE 10:e0128376. https://doi.org/10.1371/journal.pone.0128376.

Waldbusser, G. G., and J. E. Salisbury. 2014. Ocean acidification in the coastal zone from an organism's perspective: multiple system parameters, frequency domains, and habitats. Annual Review of Marine Science 6:221-247. https://doi.org/10.1146/annurev-marine-121211-172238.

Wang, N., Y.-H. Lee, and J. Lee. 2008. Recombinant perlucin nucleates the growth of calcium carbonate crystals: molecular cloning and characterization of perlucin from disk abalone, *Haliotis discus discus*. Comparative Biochemistry and Physiology B 149:354-361. https://doi.org/10.1016/j.cbpb.2007.10.007.

Wang, Q., R. Cao, X. Ning, L. You, C. Mu, C. Wang, L. Wei, M. Cong, H. Wu, and J. Zhao. 2016. Effects of ocean acidification on immune responses of the pacific oyster *Crassostrea gigas*. Fish and Shellfish Immunology 49:24-33. https://doi.org/10.1016/j.fsi.2015.12.025.

Whitehead, A., and D. L. Crawford. 2006. Variation within and among species in gene expression: raw material for evolution. Molecular Ecology 15:1197-1211. https://doi.org/10.1111/j.1365-294X.2006.02868.x.

Whitehead, A., J. L. Roach, S. Zhang, and F. Galvez. 2011. Genomic mechanisms of evolved physiological plasticity in killifish distributed along an environmental salinity gradient. Proceedings of the National Academy of Sciences of the USA 108:6193-6198. https://doi.org/10.1073/pnas.1017542108.

Wieczynski, D. J., P. E. Turner, and D. A. Vasseur. 2018. Temporally autocorrelated environmental fluctuations inhibit the
evolution of stress tolerance. American Naturalist 191:195-207.
https://doi.org/10.1086/697200.

Wigley, T. M. L., R. L. Smith, and B. D. Santer. 1998. Anthropogenic influence on the autocorrelation structure of hemispheric-mean temperatures. Science 282:1676-1679. https://doi.org/10.1126/science.282.5394.1676.

Wootton, J. T., and C. A. Catherine. 2012. Carbon system measurements and potential climatic drivers at a site of rapidly declining ocean pH. PLoS ONE 7:e53396. https://doi.org/10.1371/journal.pone.0053396.

Wootton, J. T., C. A. Pfister, and J. D. Forester. 2008. Dynamic patterns and ecological impacts of declining ocean pH in a high-resolution multi-year dataset. Proceedings of the National

Academy of Sciences of the USA 105:18848-18853. https://doi.org/10.1073/pnas.0810079105.

## References Cited Only in the Online Enhancements

Bresnahan, P. J., T. R. Martz, Y. Takeshita, K. S. Johnson, and M. LaShomb. 2014. Best practices for autonomous measurement of seawater pH with the Honeywell Durafet. Methods in Oceanography 9:44-60. https://doi.org/10.1016/j.mio.2014.08.003.

Associate Editor: Jill T. Anderson

Editor: Jennifer A. Lau

<figure><img src="image_7.png" /></figure>

"But the common black mussel, *Mitylus edulis*, and its despised neighbor, the brown horse mussel, *Modiola plicatula*, who ever saw them walk? Propulsion is not always walking. The scallop with its large adductor muscle, by snapping together its light valves, thus forcibly ejecting the water within against the water without, flits through, and sometimes even skips upon its native element, like an aquatic butterfly. But no pedestrian does so in all Mollusca-dom. Why then should not these pedate bivalves, the mussels, walk as others of their own people do? 'For want of brains!' says one. You are mistaken, sir. They have brains, the right kind too, and in the right place,—a real pedal nerve-mass, or ganglion; a little bilobed brain at the very base of the 'understanding' itself, that is, exactly under the foot, as was fabled of a very agile dancer, that his brains were in his heels." From "Mussel Climbing" by S. Lockwood (*The American Naturalist*, 1870, 4:331-336).