Ocean acidification stress index for shellfish (OASIS): Linking Pacific oyster larval survival and exposure to variable carbonate chemistry regimes

Ocean acidification (OA) results from the relatively rapid increase of partial pressure of dissolved CO₂ (PCO₂) through anthropogenic activities (Caldeira and Wickett, 2003; Orr et al., 2005; Doney et al., 2009) decoupled from concurrent changes in total alkalinity or buffering capacity that tend to happen on longer geologic time scales (Honisch et al., 2012; Zeebe, 2012). The current rate of release of anthropogenic CO₂ into the atmosphere is unparalleled in the last 66 million years (Zeebe et al., 2016), while the weathering of carbonate and silicate rocks to buffer the increased PCO₂ operates over time scales of 10,000–100,000 years (Zeebe, 2012). This temporal lag between the main global scale drivers of ocean carbonate chemistry leads to reduced pH and seawater saturation state with respect to calcium carbonate (Ω), shifting the oceans towards more corrosive conditions (Feely, 2004; Honisch et al., 2012).

Adverse biological effects associated with OA have been described across multiple stages of complex bivalve life cycles. These stages include gametogenesis (Parker et al., 2009; Barros et al., 2013), larvae (Kurihara et al., 2007; Talmage and Gobler, 2010; Barton et al., 2012; Waldbusser et al., 2015a, 2015b), juveniles (Waldbusser et al., 2011; Hettinger et al., 2012; Miller and Waldbusser, 2016) and adult calcification (Gazeau et al., 2007; Amaral et al., 2012). Emerging evidence, however, suggests an exacerbated sensitivity to OA stress during larval stages and, in particular, during the first 48 hours of bivalve larval development (Barton et al., 2012; Waldbusser et al., 2013; 2015a, 2015b; White et al., 2013, 2014). This heightened sensitivity to environmental stress during embryogenesis and early larval shell development is not restricted to acidification, but also to pollution exposure and physiological stress due to sub-optimal temperature and salinity conditions (His et al., 1999).

Bivalves are found in coastal environments and estuaries where myriad physical and biological processes result in high temporal variability in carbonate chemistry at local spatial domains (Hales et al., 2005; Feely et al., 2010; Waldbusser and Salisbury, 2014). The co-occurring modes of temporal variability in these systems include hourly, daily, weekly and seasonal time scales, resulting in complex patterns of carbonate chemistry (Harris et al., 2013; Waldbusser and Salisbury, 2014; Kapsenberg et al., 2015; Hales et al., 2017). Furthermore, the shifting baseline of atmospheric carbon dioxide concentrations (i.e., OA) will result in periods of less favorable chemistry conditions becoming more frequent and of greater magnitude (Harris et al., 2013; Hauri et al., 2013; Evans et al., 2015; Hales et al., 2017).

The physiological effects of high-frequency variability in carbonate chemistry are difficult to study and currently poorly constrained. Consequently, the incorporation of environmental fluctuations into experimental design and model predictions has been identified as a research priority within the OA research community (Takeshita et al., 2015; Wahl et al., 2016; Boyd et al., 2016). The majority of published work on bivalve larvae involves experimental designs that utilize conditions of fixed levels of high PCO₂, low pH and/or low saturation states that are maintained throughout the exposure period (Kurihara et al., 2007; Gazeau et al., 2011; Timmins-Schiffman et al., 2013; Waldbusser et al., 2015a, 2015b, 2016). In addition, the values chosen for experimental treatments are frequently based on global averages instead of local exposure regimes (McElhany and Busch, 2013), thus often not representative of the range of conditions experienced by organisms in their microhabitats (Helmuth et al., 2010; Helmuth et al., 2014).

Several recent studies have utilized different approaches to integrating the environmental variability of carbonate chemistry into manipulation experiments, aptly summarized by Boyd et al. (2016). Only a very small subset of the reviewed literature, however, concern experimental work performed on the larval stages of marine invertebrates, and currently results have been mixed in terms of whether variable or constant conditions are more favorable to organisms. In a study of two species of mussel larvae, Frieder et al. (2014) found mitigating effects of variable conditions on early development, but mixed effects on growth after eight days. Clark and Gobler (2016) also found both positive and negative effects on three different species of bivalve larvae across different physiological measures. Eriander et al. (2016), however, did not find that variable carbonate chemistry, relative to fixed, affected mean growth of barnacle larvae, although it did elicit a 20-fold increase in variance across the individuals. Overall, interpreting such complex and divergent responses could benefit from the development of conceptual and quantitative frameworks to better assess and predict the effects of dynamic acidification conditions on organisms. Understanding larval bivalve responses to exposure to variable carbonate chemistry regimes (in situ conditions) is critically dependent on realistic quantification of physiological stress on high frequency time scales and the potential for recovery after transient exposure to harmful conditions.

The effects of environmental variability on physiological stress have been noted previously (Helmuth and Hofmann, 2001), but only recently have been highlighted in OA research frameworks (Waldbusser and Salisbury, 2014) and in the limited experimental work noted above. In a seminal paper, Helmuth et al. (2010) addressed the need to focus on “organismal climatology”, or the translation of environmental conditions into realistic and relevant metrics for organisms. The exposure history of organisms that inhabit dynamic environments is driven by the “filtering” of the environmental signal or its translation into what an organism eventually experiences (Helmuth, 2009; Helmuth et al., 2010; Waldbusser and Salisbury, 2014) based on its morphological and behavioral characteristics, thus creating each organism’s fundamental physiological niche (Kearney and Porter, 2004, 2009). Organismal climatology is therefore fundamental to the development of realistic predictive models of population dynamics and species distributions (Seabra et al., 2015; Montalto et al., 2016). A wide variety of methods have been proposed, from correlative, descriptive climate-envelope models (Hijmans and Graham, 2006) to highly sophisticated multi-variable mechanistic models (Buckley et al., 2011; Montalto et al., 2016). Nevertheless, translating laboratory experimental results into mechanistic physiological models requires an abundance of physiological and environmental data, and the results are often restrictive in their application (Buckley et al., 2011). One mechanistic framework widely used across applications and relatively simple to implement is the degree-day (DD) model.

A DD model is widely used to translate the variable thermal exposure of an organism to a metric of biological performance (Bonhomme, 2000). DD models are considered parsimonious, here defined as simple parameterizations with good explanatory power, and have been used to predict crop yield as growing degree-days (GDD), to model plant development and phenology (Moore and Remais, 2014), to estimate fish growth (Neuheimer and Taggart, 2007; Chezik et al., 2014), and to forecast patterns of coral bleaching (Liu et al., 2003; Eakin et al., 2010; Kayanne, 2017). Although there are varying levels of complexity within DD models, all consider the organismal response as the result of the composite of temperature relative to a defined physiological threshold and time. For example, in DD models that predict growth, thermal exposure accumulates when temperature exceeds the minimum threshold for growth but does not surpass a second threshold for impaired physiological activity (Moore and Remais, 2014). Therefore, degree-days are referred to as “thermal time” or “thermal organismal history”.

We propose a reinterpretation of the DD framework to account for accumulated acidification stress on larval bivalves, defining a new metric: ocean acidification stress index for shellfish (OASIS). Akin to a DD model, OASIS integrates an instantaneous acidification stress effect over the larval period, is based on OA stress thresholds identified in laboratory experiments, and includes an additional parameterization that accounts for heightened sensitivity observed in early larval stages. In analogy to DD, OASIS can be described as “acidification time” or “acidification organismal history”.

We assessed the performance of OASIS by following four cohorts of Pacific oyster (Crassostrea gigas) larvae subject to a range of carbonate chemistry conditions at a commercial hatchery, Whiskey Creek Hatchery (WCH), Netarts Bay, OR, where ocean acidification has had major impacts on hatchery production of oyster seed (Barton et al., 2012; 2015). The objectives of our study were to: 1) determine differences in acute vs. chronic exposure to acidification stress; 2) document the integrated effects of exposure to buffered larval culture water vs. naturally variable carbonate chemistry on survival in four Pacific oyster larval cohorts; and 3) develop and validate the OASIS index for survival of Pacific oyster larvae to pre-settlement stage.

We measured survival to pre-settlement (i.e., pediveliger stage) of four cohorts of Pacific oyster larvae (C. gigas) exposed to different carbonate chemistry regimes at WCH: two cohorts in August 2014, one in May 2015, and one in June 2015. The experimental design was an incomplete cross-over in which C. gigas larvae were divided into two groups and reared for the first 48 hours post-fertilization, then split again into two more groups for a total of four treatments. Our designation of these treatments include capital “B” or “U”, correspoding to buffered (described below) or unbuffered during the first 48 hours, then lowercase “b” or “u” for the same treatment applied during the remaining larval period (see Table 1).

The Whiskey Creek Shellfish Hatchery (WCH) is located in Netarts Bay, OR (45.403°N, 123.944°W), a small tidally dominated estuary located on the north Oregon coast with generally near-ocean salinity in the spring and summer months (Engle et al., 2007). Although almost the entire bay is flushed with every tidal cycle (Whiting and McIntire, 1985), extensive seagrass beds (Zostera spp.) covering approximately two-thirds of the total area of the estuary (Whiting and McIntire, 1985) and microphytobenthos on tidal flats appear to significantly modify the Bay’s carbonate chemistry (Barton et al., 2012).

WCH installed carbonate chemistry monitoring equipment in 2011, revealing environmental variability on hourly, weekly and fortnightly scales (Vance, 2012; Barton et al., 2012). In 2013, the equipment was upgraded to provide real-time, high frequency measurements of PCO₂ and total CO₂ (TCO₂) that allow complete parameterization of the carbonate system (Barton et al., 2015). Previously, in 2011, WCH installed an automated system based on pH to buffer seawater, therefore enhancing its favorability for larval shell growth. WCH “buffered” water refers to incoming seawater treated with agricultural grade soda ash (unpurified sodium carbonate, Na₂CO₃) to achieve an aragonite saturation state (Ω_Ar) between 4 and 5. WCH “unbuffered” water is untreated and reflects the natural carbonate chemistry variability of Netarts Bay which modifies Pacific Ocean waters.

We constructed a flow-through experimental system consisting of sixteen 23.5-L plastic culture chambers with individual inlets reaching the near-bottom of each chamber and an outlet at the top of the culture chamber fitted with a “banjo” type filter. The banjo filter consisted of a circular section of 25-cm diameter equipped with nylon mesh screen (Nitex) to allow the flow of seawater and algae but to retain the embryos and larvae. Banjos with larger Nitex were exchanged during the experiment as larvae grew to allow better flushing and minimize clogging. The mesh size was 20 µm the first week and 60 µm from day 7 post-fertilization until the end of the experiment. The locations of the inlet and outlet were designed to maximize mixing and reduce the potential for stratification.

Experimental chambers were maintained at 25°C with four water baths (four chambers per water bath) and individual in-line heaters. The temperature of culture chambers was monitored at least twice a day. Chambers were cleaned every two or three days, coinciding with larval sampling.

The flow of seawater was individually regulated for each culture chamber with rotameter type flow meters and valves maintaining flow between 20 and 23 L hr^–1 per chamber. We estimated that the residence time of each culture chamber was less than 1.5 hours, except for a few transient occurrences of reduced or stopped flow due to routine hatchery maintenance operations (interruptions less than 2 hours long) and one instance of local disruption in electrical service (less than 8 hours).

All culture chambers were fed continuously to maintain a total approximate concentration of 25,000–30,000 cells mL^–1 of a mixture of different algal species following the practices of the hatchery (Isochrysis sp., both Caribbean and Tahitian strains, Nannochloropsis sp., Chaetoceros gracilis, and C. calcitrans).

We carried out four independent experiments on 5 August 2014, 12 August 2014, 28 May 2015 and 5 June 2015, providing four independent cohorts. For each experiment, we partnered with the hatchery managers to collect a portion of fertilized eggs from larger, commercial spawns. For each spawn, between 2 and 10 Molluscan Broodstock Program (MBP; Langdon et al., 2003) females and 1–3 MBP males were stripped of their gametes, which were then pooled and fertilized to produce one cohort.

After ~1.5 hours post-fertilization, we checked the fertilized eggs for cell division and presence of polar bodies to ensure successful fertilization. Subsequently, we divided the embryos into two groups: “buffered” (B) and “unbuffered” (U) and transferred them into eight culture chambers, four replicates per treatment, at a stocking density of 7 or14 embryos mL^–1 (Cohorts 1–3: 7 embryos mL^–1 or 165,000 embryos chamber^–1; Cohort 4: 14 embryos mL¹ or 335,000 embryos chamber^–1). Cohort 4 was stocked at a higher density to ensure sufficient survival, as this cohort coincided with background high early larval mortality experienced at the hatchery. We allowed the embryos to develop in static conditions (zero flow) for 40 to 45 hours until the Prodissoconch I shell had been fully formed and larvae had reached the D-stage. We then further sub-divided the B and U treatments into the four treatment groups (Bb, Bu, Ub, Uu) resulting in two replicate chambers per treatment. Organismal density was not adjusted for Cohorts 1–3 during the length of the test, but larval density was adjusted for Cohort 4 from 14 embryos mL^–1 to 7 larvae mL^–1 at the point (48 hours post-fertilization) when larvae were transferred to one of their four respective treatments to match the organismal density of the other cohorts during the remaining larval period. We terminated the experiments when approximately 50% of the larvae in the Bb treatment appeared to have reached the pediveliger stage and showed signs of competency to settle (i.e., presence of “eye spots” or protrusion of larval “foot”).

Larval abundance was measured at the end of the static culture period and at every chamber cleaning to measure survival. Larvae were collected on a 37-µm screen and concentrated in a known volume of seawater. At least three aliquots of equal volume were then collected and counted at 20× magnification under light microscopy. Sample volumes were adjusted to target 20–100 larvae. Counts were repeated if replicate counts differed by more than 20% among aliquots. The arithmetic mean of the counts was computed and larval survival was reported as the percentage of the total number of embryos stocked in each culture chamber. Cohort 4 survival post-48 hours was computed with respect to the corrected density of larvae.

Untreated seawater pumped straight from Netarts Bay was characterized continuously by a combined PCO₂/TCO₂ dual analyzer based on Bandstra et al. (2006) coupled with simultaneous measurements of temperature and salinity (SBE 45 MicroTSG, Sea-Bird Scientific) as described in Vance (2012). All measurements were collected at high-frequency temporal resolution (1 Hz). Independently calibrated internal gas and liquid standards ensured linearity and accuracy for both parameters.

Discrete samples for carbonate chemistry were collected from one randomly selected culture chamber of each of the two seawater sources (buffered and unbuffered) once or twice a day on larval sampling and non-sampling days, respectively. Samples were preserved in 350-mL amber glass bottles with polyurethane-lined crimp-sealed metal caps and preserved with 50 µL of saturated HgCl₂ solution. PCO₂ and TCO₂ were measured following procedures by Bandstra et al. (2006) refined for discrete samples as in Hales et al. (2005) and Hales et al. (2017). The remaining carbonate chemistry parameters were computed following water dissociation constants from Millero (1995), acid equilibration constants from Millero (2010), Dickson (1990) constants for boric acid, and calcite and aragonite solubility constants from Mucci (1983).

We opted to base OASIS on Ω_Ar based on compelling evidence that this parameter determines growth and development sensitivity to ocean acidification in early larval stages of vulnerable bivalves at modern ambient-water values (Gazeau et al., 2011; Thomsen et al., 2015; Waldbusser et al., 2015a, 2015b). Accordingly, sensitivity thresholds based on Ω_Ar have been proposed for some organisms and taxa based on fixed exposures (Barton et al., 2012; Waldbusser et al., 2015a) and used in estimations of historical environmental favorability for bivalve populations (Hales et al., 2017), biogeochemical model projections (Gruber et al., 2012), and shellfish risk assessments (Ekstrom et al., 2015).

OASIS provides a measure of accumulated stress by integrating only the area of the Ω_Ar time series that falls below the designated threshold over the entire larval period, from fertilization to metamorphosis. Mathematically, OASIS can be described as follows:

where Ω_thrshd is the selected acidification stress threshold, Ω_t is the instantaneous value, dt is the time differential (minutes), D_pft is day(s) post-fertilization (spawning day = 0), and if (Ω_thrshd – Ω_t) > 0 is the truncating condition of the integral, allowing the integral to operate only when the conditions are below the set threshold. The units of OASIS are Ω min day^–1. It follows, therefore, that OASIS will have different total integration times depending on developmental rates. Longer periods from fertilization to pre-settlement will result in increased integration time, but only in increased OASIS values if Ω_thrshd > Ω_t.

The most important parameterization of our model is the selection of a physiological relevant threshold. We chose a value of Ω_Ar = 1.5 as a conservative stress threshold for OASIS based on recently published work (Barton et al., 2012; Gruber et al., 2012; Ekstrom et al., 2015) that identified this value as a minimum threshold for early development and commercial production-to-settlement competency for C. gigas. The second parameterization within OASIS incorporates the increased sensitivity observed in the early larval stages of C. gigas by introducing a weight function into the metric, as described by the term $11+Dpft$⁠. This term describes an inverse and injective function with maximum of 1 and horizontal asymptote that approaches 0 when D_pft > 0; therefore, it works as a simple but useful weight function by placing maximum weight on the early larval stages (D_pft = 0 corresponds to a weight score of 1; D_pft = 2, weight score of 0.33; D_pft = 7, weight score of 0.125).

We computed OASIS for the cohorts described above by combining results of discrete carbonate chemistry sampling for the static portion of the experiment (i.e., embryogenesis to Prodissoconch I shell) with continuous, high-frequency high-resolution data collected by the PCO₂/TCO₂ analyzer for the remainder of the experiment. As the embryos were raised in static conditions for the first 48 hours, we computed the OASIS value for the unbuffered treatment by fitting a linear equation between the initial and final Ω_Ar calculated from the results of discrete carbonate chemistry sampling of the embryo culture chamber. We acknowledge that Ω_Ar likely varied in a non-linear fashion during these time periods due to respiration and calcification, but decided against taking excessive discrete samples to avoid causing stress to the embryos during this critical life stage. If Ω_Ar was above the defined threshold of Ω_Ar = 1.5 in the discrete culture chamber, then the OASIS value was 0 during this period. After 48 hours, and thus during the period of flow-through culture, the continuous record was used to calculate OASIS in the unbuffered treatments. OASIS was zero in all buffered treatments for all periods because Ω_Ar was always greater than 1.5.

Our experimental design was geared towards an analysis of variance (ANOVA) of larval survival as a function of treatment; however, the survival data did not meet ANOVA assumptions of normality and homoscedasticity even after multiple data transformations were tested. We therefore explored the use of generalized linear models.

Fractional linear regression (Papke and Wooldridge, 1996), a type of generalized linear model, is suitable for fractional responses contained between 0 and 1. It approximates the sample distribution as logistic, allows for overdispersion of the data, and uses iterative weighted least squares to fit a linear model. Also, to address the heteroscedasticity of the data, we computed heteroscedasticity-consistent (HC) standard errors, also known as robust or White-Huber standard errors, and compared them to the regular standard error (SE). Specifically, we used HC3, a version particularly well suited for small sample sizes (Long and Ervin, 2000). Differences between regular and robust standard errors were deemed insignificant, as they did not change statistical significance in any of the parameters, so we decided to report regular SE.

We analyzed survival data at 48 hours post-fertilization and pre-settlement independently, given the autocorrelation of observations over time. All cohorts were grouped for each analysis, using cohort as an independent factor. Therefore cohort had four levels and treatment had two levels (buffered, B; unbuffered, U) for the analysis at 48 hours post-fertilization and four levels (buffered-buffered, Bb; buffered-unbuffered, Bu; unbuffered-buffered, Ub; unbuffered-unbuffered, Uu) for the analysis at pre-settlement. In all cases, we also tested the interaction of cohort and treatment.

We used standard linear regression to fit survival data to OASIS and we evaluated the statistical significance of the correlation with Pearson’s coefficient. All statistical analyses were performed in Rstudio (version 1.0.136).

The buffered seawater treatment consistently had greater Ω_Ar values than the unbuffered one for all cohorts (Tables 2 and 3; Figures 1 and 2). Mean conditions for unbuffered seawater carbonate chemistry during the static culture period (first 48 hours post-fertilization) were below our stress threshold only for Cohort 2 (Ω_Ar = 1.34), close to threshold for both Cohorts 1 and 4 (Ω_{A r} = 2.25 and Ω_Ar = 1.87, respectively), and well above the threshold for Cohort 3 (mean Ω_Ar = 3.22) (Table 3). For the remainder of the larval period, the unbuffered carbonate chemistry conditions experienced by Cohorts 3 and 4 were less acidified than for Cohorts 1 and 2, as evidenced by the reduced overall mean values of PCO₂ (1288 µatm and 882 µatm for Cohorts 3 and 4, compared to 1445 µatm and 1520 µatm for Cohorts 1 and 2, respectively) and, conversely, overall greater values for pH (7.76 and 7.79 for Cohorts 3 and 4, versus 7.60 for Cohort 1 and 7.54 for Cohort 2) and Ω_Ar (2.14 and 2.13 for Cohorts 3 and 4, compared to 1.38 for Cohort 1 and 1.18 for Cohort 2). The mean unbuffered Ω_Ar was slightly lower than the defined threshold for physiological stress for Cohorts 1 and 2 (Ω_Ar < 1.5), but slightly higher for Cohorts 3 and 4 (Ω_Ar > 2.0). Buffered treatments were always above the physiological threshold (Ω_Ar > 4.0) (Table 3). The temperature remained reasonably stable around 24.3°C (± 0.8°C standard deviation, STD) due to the controlled heating of water in the hatchery, and salinity was generally high (33.0 ± 0.9 STD).

Our culture chambers reflected the dynamics of the source water to the hatchery as evidenced by the general agreement of the discrete chamber samples with the continuous high-resolution carbonate chemistry record (see Figures 2 and 3). Small discrepancies, more evident in PCO₂ and pH than in Ω_Ar, were likely the result of a combination of processes including the heating of the water to 24.4°C, the interplay of larval calcification and net respiration, and the addition of concentrated CO₂-enriched algal cultures.

Larval survival at 48 hours post-fertilization and pre-settlement are shown in Figures 3 and 4. Analysis of all cohorts together yielded a highly significant effect of treatment on survival to 48 hours post-fertilization (F_1,28 = 51.51, p = 2.04 × 10^–7), while cohort (F_3,28 = 7.05, p = 0.0015) and the interaction between cohort and treatment (F_3,28 = 9.90, p = 0.00020) were also significant, limiting overall inferential power. Survival to 48 hours post-fertilization varied among treatments and across cohorts (Figure 3, Table 3). Cohort 1 showed a strong survival response to buffering (B = 88.6 ± 7.8% SE versus U = 7.6 ± 8.6% SE), while Cohorts 2 and 4 showed possible trends towards increased survival (B = 53.7 ± 4.4% SE versus U = 38.5 ± 13.5% SE for Cohort 2; B = 96.0% ± 6.9% SE and U = 85.5% ± 15.3% SE for Cohort 4). In contrast, Cohort 3 displayed a possible reverse trend, though the mean survivals were similar (B = 79.3% ± 6.1% SE versus U = 81.5% ± 7.5% SE; Figure 3, Table 3).

Three of the four cohorts, however, showed strong buffering effects in pre-settlement survival (Figure 4). We decided to exclude Cohort 4 from the analysis, after suspicions that other factors beyond our treatments had affected the survival responses (see Discussion). After excluding Cohort 4, survival was influenced by both treatment (F_3,23 = 3.39, p = 0.041) and cohort (F_2,23 = 4.91, p = 0.020), but not their interaction. However, none of the post-hoc Tukey pair-wise comparisons were significant, due to the overall high variance; therefore, any observed differences have to be interpreted with caution.

As expected for all cohorts, we found that Cohorts 1 and 2 showed higher survivorship for buffered larvae, compared to those exposed to naturally unbuffered conditions. For Cohort 1 only those larvae buffered during the early larval stage (i.e., first 48 hours post-fertilization) survived to pre-settlement stage (5.4% ± 5.4% SE and 1.6% ± 0.6% SE, for Bb and Bu, respectively, versus 0.0% ± 0.0% SE for all unbuffered). In contrast, Cohort 2 survival was highest for those larvae buffered during the whole experiment (12.9% ± 1.9% SE) and lowest for those unbuffered (5.5% ± 5.5% SE), but was similar for those buffered only for the first 48 hours (8.1% ± 3.1% SE) and those buffered only after this period (8.1% ± 2.2% SE). For Cohort 3, the pattern observed for Cohort 1 was repeated: larvae buffered during early stages outperformed those that were not buffered, although Cohort 3 seemed to benefit from exposure to natural variability after initial buffering (11.1% ± 7.8% SE for Bb; 18.6% ± 2.0% SE for Bu; 1.2% ± 0.9% SE for Ub; and 1.9% ± 1.7% SE for Uu). Finally, Cohort 4 showed an unexpected pattern of survival, in which buffering the water seemed to have a detrimental effect, and survival was highest in larvae not exposed to buffered water at any point during the experiment (3.4% ± 0.3% SE for Bb; 2.6% ± 2.6% SE for Bu; 1.9% ± 1.4% SE for Ub and 12.0% ± 2.6% SE for Uu).

Competency to settle in the fully buffered (Bb) treatment was reached at different days post-fertilization for each cohort, resulting in varying duration of the experiments. Cohort 1 did not reach pediveliger stage until 22 days post-fertilization, whereas Cohort 2 reached competency to settlement at 16 days, and Cohorts 3 and 4 at 17 days (Table 3).

The Ω_Ar time series for Cohorts 1, 2 and 3 and the associated computed OASIS values are shown in Figure 5. We excluded Cohort 4 from OASIS calculations given the odd pattern of survival (see Discussion).

The unbuffered Ω_Ar during the first 48 hours was above the stress threshold for Cohorts 1 (mean Ω_Ar = 2.25) and 3 (mean Ω_Ar = 3.22), so the accumulated OASIS during this stage is 0 for all treatments of these two cohorts (Figure 5A, C). For Cohort 2, we assumed that Ω_Ar increased linearly from Ω_Ar = 1.00 to Ω_Ar = 1.67 during the first 48 hours (Figure 5B). The cumulative OASIS value for Bb for all three of these cohorts was 0, as Ω_Ar never reached values below the stress threshold.

Given that OASIS was 0 for the first 48 hours for Cohorts 1 and 3, OASIS values for Bu and Uu were identical: 741 Ω min day^–1 for both values for Cohort 1 and 253 Ω min day^–1 for Cohort 3. In the case of Cohort 2, the OASIS values for Bu, Ub, and Uu were 499, 557, and 1056 Ω min day^–1, respectively. The integration period was longer for Cohort 1 (22 days) compared to Cohorts 2 and 3 (16–17 days) as a consequence of its slower development.

OASIS predicted survival for 2 of the 3 cohorts tested. The regression fit between survival and OASIS for Cohorts 1, 2, and 3 is shown in Figure 6. We excluded from the computation two treatments from Cohort 1 (Ub and Uu) due to concern about the validity of a cumulative stress index to assess pre-settlement survival when no survival to the end of the larval stage was achieved. Survival reached 0% at 2 and 5 days post-fertilization for Cohort 1 Uu and Bu, respectively. As two data points have been excluded in the case of Cohort 1, the relationship is a perfect fit (%Survival = –0.005[OASIS] + 5.352). We recognize that fitting a linear regression to the remaining two points is problematic, and only provide the values to compare the slopes of the two cohorts. We found a significant linear regression for Cohort 2 (%Survival = –0.007[OASIS] + 12.346; R² = 0.957, p < 0.05), but no correlation for Cohort 3 (R² < 0.1).

This study was designed to address three main objectives related to the effects of variable carbonate chemistry on larval oyster survival. First, we assessed the effect of acute exposure (first 48 hours post-fertilization) to naturally occurring carbonate chemistry on larval survival and compared it to larvae cultured in buffered seawater. Second, we determined chronic effects of buffering as an OA mitigation strategy for the entire larval stage. Third, we tested the performance of OASIS as a predictor of survival to settlement competency on larval stages exposed to different carbonate chemistry regimes.

Inquiry into the effects of buffering as a mitigation strategy for acute exposure to acidification stress (first 48 hours post-fertilization) revealed that high Ω_Ar during embryogenesis dramatically increased survival to 48 hours post-fertilization in Cohort 1, with either a possible trend towards increased survival or else no effect detected in the other cohorts (Figure 3). These results appear generally consistent with earlier work (Gazeau et al., 2011; Thomsen et al., 2015; Waldbusser et al., 2015a, 2015b), showing reduced acute sensitivity in development and growth to ambient-water Ω_Ar values exceeding 1.5. However, comparisons with previously published work are difficult, as survival is a different metric than development and growth. Further complicating comparisons, some experimental studies have been conducted at temperatures that span the range of sub-optimal conditions for larval growth and survival (Bochenek et al., 2001).

Cohorts 1 and 4 had very similar chemistry conditions of about Ω_Ar = 2.0 during the first 48 hours yet showed remarkably dissimilar responses. Secondary effects from other carbonate chemistry parameters are unlikely to have been responsible for these divergent results, given that PCO₂ was lower for Cohort 1 than for Cohort 4 (Table 3). Although no definitive studies on survival are available, robust experimental data have demonstrated the sensitivity of development and growth to Ω_Ar, but not to PCO₂ and pH at the embryological stage (Waldbusser, 2015a). We propose, therefore, that the strong response of Cohort 1 to buffering might be related to additional environmental stressors or to broodstock quality that resulted in an overall poorer condition. Cohort 1 larvae might have had greater sensitivity to acute exposure at higher Ω_Ar if their overall fitness had already been compromised. The overall longer larval period, suggesting developmental delays (as in Timmins-Schiffman et al., 2013), might also be evidence for reduced fitness.

Buffering seawater during embryogenesis resulted in significantly greater pre-settlement survival in three of the four cohorts (Figure 4). Buffering only after this early stage seemed to improve survival to pre-settlement in only one cohort when compared to survival in continuously unbuffered seawater (Figure 4B), highlighting again the heightened sensitivity to early exposure and the importance of delayed effects.

We hypothesize that the unbuffered Ω_Ar conditions experienced by Cohort 3 during embryogenesis (Ω_Ar > 3.0) were not sufficiently unfavorable to elicit an acute mortality response, but may have resulted in sub-lethal effects. Such delayed effects are evidenced by an order of magnitude greater survival to pre-settlement for larvae that were buffered during this early stage compared to those that experienced unbuffered conditions (Figure 4C). These results also confirm observations made by hatchery managers and in other studies of carryover effects that the impacts of early exposure to acidification stress often manifest themselves later in the larval cycle (Barton et al., 2012, 2015; Hettinger et al., 2012).

The overall beneficial effect of buffering agrees with observations made by multiple hatchery operators on the US West Coast (Barton et al., 2015). Better bivalve larval performance at elevated Ω_Ar (Ω_Ar > 1.5) during embryological and early larval development has also been shown in multiple laboratory experiments (Gazeau et al., 2011; Timmins-Schiffman et al., 2013; White et al., 2013, 2014; Thomsen et al., 2015; Waldbusser et al., 2015a, 2015b). In this study, we exposed some of the larvae to acute stress followed by chronic stress (Uu). The overall response of larvae to chronic stress in the Uu treatments may represent a conservative estimate of the effects of chronic exposure, as we were evaluating only the survivors of acute stress that might have been selected for increased tolerance to stress.

The overall survival-to-settlement competency size in this study fell below commercial hatchery targets (30–60%; A Barton, personal communication; Breese and Malouf, 1975) for all cohorts, but was within the range of other published studies for oyster larvae (Talmage and Gobler, 2012; Eierman and Hare, 2013; Ko et al., 2013). To our knowledge, however, laboratory studies that have measured long-term larval survival under acidification have not captured embryological exposure to acidified conditions, only the exposure of larvae at least 24 hours old (Talmage and Gobler, 2010; Ko et al., 2013, 2014), potentially underestimating deleterious effects. While our results generally show better survival for larvae raised in buffered water, differences in survival are complicated to assess statistically given the limited number of replicates due to logistical constraints, including limited laboratory space that prevented the installation of more culture chambers, lack of feasibility of sampling efforts, and the prioritization of repetition of the experiment (i.e., more cohorts) over increased sample replication. Our unbuffered treatment was also a function of natural variability, making the magnitude of the acute stress exposure variable by design.

Cohort 1 performed better under buffered conditions, but showed lower overall survival across treatments, both at 48 hours post-fertilization and pre-settlement, compared to the other three treatments (Figures 3 and 4). This reduction in overall fitness suggests possible parental effects, including lower maternal lipid contribution to eggs (Gallager and Mann, 1986) or an unfavorable condition of the broodstock (Utting and Millican, 1997). The overall poorer performance of Cohort 1 compared to the other cohorts is also evident by the notably longer time to reach settlement competency (22 days, compared to 16–17 days for the other cohorts). Although arrested or slow larval development might sometimes be an effective strategy to survive stress (Moran and Manahan, 2004), in general a longer pelagic larval stage increases the probability of loss by predation (Underwood and Fairweather, 1989) or dispersal away from favorable settlement substrates and habitats (Strathmann, 1985). Therefore, developmental larval delays can have significant detrimental effects on survival to post-metamorphosis and juvenile growth (Pechenik, 1990; 2006) and, ultimately, affect population dynamics (Pechenik, 1999).

Cohort 4 responded differently in pre-settlement survival with respect to buffering than the other three cohorts, resulting in highest survival in the continuously unbuffered treatment (Figure 4D). In 2015, an unprecedented harmful algal bloom (HAB) of diatoms, Pseudo-nitzschia spp., producers of neurotoxin domoic acid (DA), was detected along the West Coast of the US (Cavole et al., 2016; Du et al., 2016; McCabe et al., 2016; McKibben et al., 2017). One particularly severe peak coincided with the culture of Cohort 4 in early June 2015 (Du et al., 2016; Figure 4 in McCabe et al., 2016), coinciding with low embryo survival and arrested growth of later larval stages at WCH. The concentrations of particulate DA recorded during this HAB (~20,000 ng mL^–1 in California (McCabe et al., 2016)) were 2–3 orders of magnitude higher than levels reported to cause impaired bivalve larval development, growth, survival and immune system functioning (De Rijcke et al., 2015; Liu et al., 2007). Various non-mutually exclusive mechanisms to trigger DA production and toxicity have been proposed (reviewed by Lelong et al., 2012; Lewitus et al., 2012), including elevated production of DA with increasing pH (pH from 7.9 to 8.9; Lundholm et al., 2004; Trimborn et al., 2008; but see Sun et al., 2011) and greater availability of Fe in Fe-depleted acclimated cells (reviewed by Trainer et al., 2012). We hypothesize that the negative response to buffering observed in Cohort 4 could be associated with an increased production of DA in buffered water due to higher pH (mean pH_unbuffered = 7.79 versus mean pH_buffered = 8.15; Table 2), and with increased levels of Fe associated with the WCH-buffered water (Waldbusser GG, limited unpublished data), possibly due to impurities in the agricultural grade Na₂CO₃ used as buffering reagent. These observations suggest the need for better evaluation of the efficiency of routine buffering practices at commercial larval hatcheries during HABs. Importantly, if the increased pH due to buffering were to result in greater DA levels during HABS, this outcome would suggest a possible narrowing of effective ocean acidification mitigation strategies during HABs, which are predicted to increase in frequency and duration (McKibben et al., 2017). At least for early larval stages, current studies suggest that the development of buffering strategies that increase Ω_Ar while limiting changes in pH (Waldbusser et al., 2015a; 2015b; 2016) would be useful to mitigate OA effects on larval development without stimulating increases in DA production.

The ultimate objective of this study, however, was to present a framework to translate temporal variability of environmental carbonate chemistry into a physiological stress index. OASIS adequately predicted survival due to variable acidification exposure for two cohorts of Pacific oyster larvae in 2014 (Figure 6). Given the high Ω_Ar corresponding to exposure during embryogenesis for Cohorts 1 and 3, Bu and Uu values were equal. In the case of Cohort 2, similar values for Bu and Ub are striking and suggest that the accumulated acidification stress only for the first 48 hours was as severe as the cumulative stress from day 2 to pre-settlement stage.

Comparing the correlations between survival and OASIS for both Cohort 1 and 2 reveals that the slopes of the regressions are very similar, but the intercepts are different (Figure 6). Although assessing the intercept of Cohort 1 requires caution, given that only two points were included in the regression, we hypothesize that the differences in intercept may be hinting at general fitness of the cohorts. Fundamentally, the intercept models survival in the absence of acidification stress (Bb). The smaller intercept for Cohort 1, which performed worse than the other three cohorts (Figures 4 and 5), might be another piece of evidence pointing towards overall lower cohort fitness possibly attributable to broodstock quality (Utting and Millican, 1997) or maternal lipid investment (Gallager and Mann, 1986). This observation suggests that besides predicting survival associated with exposure to variable acidified conditions, OASIS might be able to help identify differences in overall performance among cohorts exposed to similar variability in the natural environment, aiding in identifying particular strains and families with better resistance to acidification for genetic selection (de Melo et al., 2016).

An advantage of OASIS is that it offers a flexible parameterization that can easily incorporate multiple thresholds across larval stages, potential harmful effects associated with high rates of change in environmental conditions, and physiological recovery. For this study, we chose a conservative physiological stress threshold (Ω_Ar = 1.5) already applied to other studies (Gruber et al., 2012; Ekstrom et al., 2015; Hales et al., 2017) and derived from hatchery production observations (Barton et al., 2012) and laboratory experiments measuring early development and growth (Waldbusser et al., 2015a, 2015b). Survival, however, may be more sensitive to reduced Ω_Ar conditions, particularly in early larval stages, as evidenced by the contrasting responses observed with Ω_Ar values around 2.0.

Although we account for changes in sensitivity through a weighting function, OASIS currently assumes a constant stress threshold over the entire larval stage. A secondary threshold that would result in a certain level of cumulative sub-lethal stress when early larval exposure Ω_Ar is > 1.5 and < 4.0 might be needed, as evidenced by the reduced pre-settlement survival of Cohort 3 larvae exposed to Ω_Ar > 3.0 during the first 48 hours. More work is needed to assess if the results obtained with Cohort 3 are the result of a new secondary threshold for this larval stage. Multiple thresholds, however, could be incorporated to further refine stress exposures over the entire larval cycle, or to include the presence of multiple known or unknown stressors (for example, HABs in this study and Talmage and Gobler, 2012; hypoxia in Keppel et al., 2015, and Clark and Gobler, 2016).

Additional terms can also be added to OASIS to incorporate additional stress due to highly dynamic environmental conditions like those experienced by Cohort 3 at eight days post-fertilization when Ω_Ar decreased from 3.39 to 1.18 in less than 2 hours (Figure 5C). Some organisms have already been found to be sensitive to rapidly changing conditions in laboratory experiments. For example, Kamenos et al. (2013) demonstrated that coralline algal structure is more sensitive to a rapid change of PCO₂ conditions than to the magnitude of the perturbation.

Finally, the potential for physiological recovery during times when carbonate chemistry is more favorable has not been adequately explored in the literature yet. Currently, OASIS only integrates stressful conditions; it is, therefore, a truncated integral. Similarly, we lack enough data to parameterize the potentially additive deleterious effect of persistent harmful conditions. As more data regarding physiological effects of variable carbonate exposure are collected, recovery thresholds and parameters to account for continuous unfavorable chemistry can be implemented in OASIS. At this stage in its development, the current parsimonious description of OASIS predicted survival in two of the three cohorts tested without adding more layers of complexity. The effectiveness of OASIS for the other two cohorts appears to be limited not by the approach but by the lack of OA stress in the first 48 hours for Cohort 3 and by an unprecedented HAB that appeared to overcome any OA effects, including the possibility of increased toxicity in the buffered conditions, for Cohort 4. Therefore, OASIS represents an encouraging advance towards the task of translating laboratory results into more meaningful metrics that accurately capture accumulated physiological stress due to organismal exposure to variable acidification conditions, but more work is needed to further validate its performance and better understand the role of multi-stressors.

We believe that OASIS, a re-interpretation of DD models, is a promising step towards understanding the effects of variable carbonate chemistry experienced by bivalve larvae in the field. More data regarding experimental exposure to fluctuating treatments are needed to better understand the underlying physiological mechanisms of sensitivity to variable carbonate chemistry and to reconcile conflicting published results (Frieder et al., 2014; Keppel et al., 2015; Eriander et al., 2016; Clark and Gobler, 2016), but OASIS is a flexible framework that can continue to be refined to distill the increased complexity in responses and facilitate the prediction of outcomes (Helmuth et al., 2014; Boyd et al., 2016). Multiple strategies should be pursued simultaneously: continued monitoring efforts to characterize and analyze environmental conditions in temporal and spatial scales relevant to organisms (Helmuth et al., 2010; Dillon and Woods, 2016); development and validation of mechanistic models that incorporate physiological sensitivity into biological performance (this study; Kearney and Porter, 2009; Helmuth et al., 2014; Bozinovic and Pörtner, 2015); and integration of these models into local, regional and global climate models to assess the future of populations and their distributions (Buckley et al., 2011).

Although a growing body of experimental work has advanced our understanding of the effects of ocean acidification on marine bivalves, key vulnerable taxa, important questions remain regarding the effects of highly variable carbonate chemistry commonly found in bivalve habitats on larval stages and the effectiveness of buffering mitigation strategies in current use within commercial larval hatcheries. Our work supports the benefits of buffering seawater, particularly during embryogenesis and the formation of the prodissoconch I shell, on overall larval survival. Increasing coverage of high frequency monitoring of carbonate chemistry is demonstrating that carbonate chemistry is dynamic in coastal environments. Translating those variable data into reasonable approximations of success for sensitive species and taxa is important for predicting OA effects on populations and communities. To help translate the effects of dynamic acidification conditions on bivalve larvae survival, we propose a new quantitative framework, the ocean acidification stress index for shellfish (OASIS), based on a conservative threshold derived from past observations that accounted for differential sensitivity across the larval stage. OASIS predicted survival for two of the three larval cohorts exposed to different carbonate chemistry regimes, thus representing an important step towards distilling physiological responses to complex carbonate chemistry conditions into synthetic metrics useful to industry managers and policy-makers as predictive tools.

Biological data, OASIS computation results, and discrete carbonate chemistry data are available at the UC Press Dash online repository (https://doi.org/10.15146/R3VT1Z). Continuous carbonate chemistry monitoring data are available through the website portals for NANOOS (http://www.nanoos.org/) and GOA-ON (http://portal.goa-on.org/) or upon request to Dr. Burke Hales.

The authors thank the Whiskey Creek Hatchery personnel for providing space, broodstock, and countless support in conducting the work within the hatchery. Particularly, the authors would like to thank Sue Cudd, Mark Wiegardt, and Alan Barton from WCH for their support, encouragement and professional experience.

This work was supported by Oregon Sea Grant grant #NA10OAR4170059 awarded to GGW and BH. IG completed this work as partial fulfillment of Ph.D. degree requirements.

The authors have no competing interests to declare.

• Contributed to conception and design: IG, GGW, BH

• Contributed to acquisition of data: IG, GGW, BH

• Contributed to analysis and interpretation of data: IG, GGW, BH

• Drafted and/or revised the article: IG, GGW, BH

• Approved the submitted version for publication: IG, GGW, BH