Systematic Review and Meta-Analysis Toward Synthesis of Thresholds of Ocean Acidification Impacts on Calcifying Pteropods and Interactions With Warming

Introduction

Ocean acidification (OA) is already impacting marine habitats worldwide, with projected changes to become widespread over this century (Bakun et al., 2010). The greatest rate of change is anticipated in the upwelling and polar regions (Bopp et al., 2013; Chavez et al., 2017), where OA is often compounded by co-occurring stressors, such as hypoxia, thermal stress, and increasing pollution sources from urbanizing coastlines. Concern is growing on the possibility that combined stressors may push ecosystems beyond their tipping points and deleteriously impact ecosystem services (Groffman et al., 2006; Crain et al., 2008; Halpern et al., 2008; Barbier et al., 2011). The characterization of ecosystem vulnerability requires knowledge of OA thresholds at which lethal vs. sublethal impacts occur; however, despite rapidly growing literature documenting the biological impacts of OA, scientific consensus on these thresholds is lacking. These thresholds are urgently needed to provide consistent interpretation of monitoring data and regional oceanic model output, as well as to support marine climate change vulnerability assessments. They are also an important precursor to establishing biologically-relevant OA Water Quality Goals (WQG; Aminzadeh et al., 2016; Weisberg et al., 2016).

Pteropods are appropriate indicators for such an assessment for several reasons. First, they are ubiquitous holoplanktonic calcifiers that efficiently transfer energy from phytoplankton to higher trophic levels (Lalli and Gilmer, 1989; Hunt et al., 2008), serving as an important prey group for ecologically and economically important fishes, bird, and whale diets (Armstrong et al., 2005; Aydin et al., 2005; Karpenko et al., 2007). Second, they significantly contribute to carbonate production (Bednaršek et al., 2012a). Third, they have a well-documented and specific sensitivity to OA, ranging from evidence of exposure (e.g., shell dissolution), sublethal (without lethal outcome, e.g., growth, calcification, gut clearance) responses (Lischka et al., 2011; Lischka and Riebesell, 2012; Bednaršek et al., 2017; etc.; Figures 1, 2), occurring through a variety of impairment processes (Figure 1). Experimental results with negative pteropod responses are, in majority of studies, demonstrated against low Ω_ar, and only in few studies against pH or pCO₂. Strong co-linearity among the carbonate chemistry parameters and a subsequent lack of decoupled carbonate chemistry parameters in the experimental set-ups makes it difficult to delineate the exact mechanisms of impairments (Figure 1). In general, Ω_ar or the availability of ${\text{CO}}_3^{2-}$, is a driver behind reduced biomineralization processes and growth, CO₂ -induced hypercapnia, and low pH alters metabolisms and acid-base balance, while pH seems to be involved in the neurologically induced behavioral changes, and reproductive impairments (Manno et al., 2016; Moya et al., 2016). Overall survival is likely impacted in a cumulative way by a combination of various sublethal processes, and thus not necessarily linked to a single carbonate chemistry parameter (sensu Bednaršek et al., 2017).

Figure 1. Responses of pteropods to OA. Effects are mediated via diffusive CO₂ entry (black arrows) into body and cell compartments, resulting in a rise in intracellular pCO₂ (highlighted in red), a drop in compartmental pH (highlighted in blue), and their effects (red arrows) on various biological endpoint (red text) in tissues and cellular compartments. In addition, lower aragonite saturation state (Ω_ar) and carbonate ions (${\text{CO}}_3^{2-}$) has an effect on shell dissolution (depicted as pteropod with intact and damaged shell) and reduced calcification. Figure modified from IPCC (2014).

Figure 2. (A) Summary of experimental studies with respect to the Ω_ar conditions and (B) the duration of experimental exposure used in treatments.

For these reasons, pteropods have been proposed as OA indicators for climate change assessments, ecosystem-based management, and the establishment of WQG (Weisberg et al., 2016; Bednaršek et al., 2017). Comprehensive reviews have recently summarized knowledge of pteropod vulnerability to either OA (Manno et al., 2016) or multiple stressors (Bednaršek et al., 2016). However, despite numerous experimental and field stress-response studies, no synthesis has derived pteropod-specific OA thresholds, which are key to predictions of climate change impacts, habitat suitability, and climate assessment.

Interaction among multiple stressors, where the effect of one stressor is co-dependent on the other stressors, can lead to nonlinear outcomes, such that the combined effects can have outcomes that are more severe than expected. These combined impacts are dependent on how the stressors reinforce (i.e., synergistic) or counter (i.e., antagonistic) one another. Understanding how these responses may change when co-occurring with other climate change stressors, i.e., temperature or low dissolved oxygen (DO), is key for constraining the thresholds and their application. While meta-analyses across different functional groups mostly incorporate one stressor (Dupont et al., 2010), only a few have considered the multiple stressor effect (Kroeker et al., 2013), demonstrating that meta-analyses are an important analytical tool to expand our understanding of pteropod vulnerability to multiple stressors.

This study presents a global literature review, meta-analysis, and consensus of a panel of pteropods experts on thresholds for lethal and sublethal responses to OA stress for calcifying pteropods. Further, using meta-analysis, the magnitude of thermal stress and OA on pteropod responses was investigated for the interactive effect of two stressors potentially changing these thresholds (Figures 1, 2). The review further discusses considerations for the application of thresholds to monitor data and model output in the context of marine climate change vulnerability assessments. Metrics of exposure, such as duration (time interval below the threshold), intensity (magnitude of departure from threshold), and severity were recommended for consistent application across monitoring data and model output (Hauri et al., 2013). To demonstrate their intended use, thresholds were applied to a numerical ocean model simulation from the southern California Current System (CCS). Finally, specific recommendations are provided to increase the confidence in threshold application to demonstrate their applicability in the marine climate change vulnerability assessments.

Materials and Methods

The synthesis was conducted in three parts. The first part was a meta-analysis to determine the extent to which responses among studies conducted on congener species, differing life stages, and disparate geographies could be integrated into a common analysis. The second was a breakpoint analysis to identify thresholds of OA stress for several endpoints. The third was use of professional judgment among pteropod experts to assign an uncertainty level to each threshold based on the data quantity, data quality, and consistency in findings among studies included in the breakpoint analysis. The experts also identified the exposure duration most appropriate to thresholds for each pathway.

Meta-Analysis

The global literature search and review identified 15 pteropod studies (Table 1), based on a minimum criterion of (1) experimental studies with a minimum of three Ω_ar treatments (control, mid- and high range) or (2) field stress-response studies conducted across a range of Ω_ar values (e.g., Bednaršek et al., 2014b, 2017). Such criteria are needed because the meta-analyses results provided the base for the threshold analyses; hence, the same studies needed to be evaluated in both meta- and breakpoint-analyses, excluding the studies with only two different treatments. The selection was based on the available datasets that were published up to June 2017, with 9 different species across different areas (west coast, temporal, tropical, and polar) and 22 different pteropod endpoints recorded in the studies. The exception for this was meta-analyses with two stressors, OA and warming, in which the results were not used for threshold analyses, with all of the available data from literature until June 2017 specifically related to Limacina helicina. Studies with temperature and OA reference is included in Table 2. The endpoints were characterized into three effect categories, ranging from indications of exposure response, sublethal and lethal response across various species, life stages, and geographically different locations. A total of 2,097 datapoints were extracted from the raw experimental data provided in published literature or from authors with unpublished data. While some studies used multiple exposure metrics, every one of the studies used aragonite saturation state (Ω_ar) as the primary exposure metric; as such, all data were placed in the context of this single parameter for this study. Meta-analysis was used to (1) evaluate the robustness and signal-to-noise ratio of pteropod response measures relative to OA stress and (2) evaluate the extent to which disparate endpoints or experimental datasets can be combined given the variation in potentially confounding factors, such as intraspecific differences, the latitudinal ranges, temperature and DO gradients, life stage, duration of treatment, and nutritional/energetic status of experimental organisms. Meta-analysis of pteropod OA response measures generally followed the methods of Arnqvist and Wooster (1995), but the quantitative assessments of the effect size were derived from raw experimental data, in contrast with previous studies using a meta-analysis approach (e.g., Harvey et al., 2013; Kroeker et al., 2013; Meyer and Riebesell, 2015). First, the OA effect observed in all experiments was assessed using a unified approach. For each response measure, the experimental values (Y) were plotted against aragonite saturation state (Ω_ar). Two Ω_ar thresholds were selected by optimizing the differences between the “Effect” (low Ω_ar) and the “Control” (high Ω_ar) groups of observations. This optimization was achieved by calculating t-statistics of the differences between the “Effect” and “Control” groups across the range of possible thresholds (Supplementary Figure 3). Data points falling between the two Ω_ar thresholds were excluded from subsequent calculations. For the “Control” group, the averaged values of the experimental measurement (Ycontrol) and saturation state (Ω_ar) were calculated. Second, all data points in the “Effect” group Yeffect[i] were transformed into a “common currency,” called “Logarithmic response ratios” (LnRR) calculated as LnRR[i] = log (Yeffect[i] / Ycontrol; Hedges et al., 1999), and the mean and standard deviation of saturation state of the “Effect” group (Ω_ar) were calculated. Significance of the OA effect for each measurement type was assessed using a Student t-statistic (difference from zero at 95% confidence interval) and the differences between the OA effects on various measurement types were demonstrated using forest plots (Lewis and Clarke, 2001). The significance level of the total OA effect, summarizing all response measures, was calculated using the R package “metaphor” (Viechtbauer, 2010). A random-effects model was used, based on the assumption that the OA effect varies between studies. The total effect was assessed as an averaged LnRR with a 95% confidence interval and the null hypothesis of the absence of the OA effect was tested. The total heterogeneity between the response measures was assessed by the coefficients T2 (total amount of heterogeneity) and I2 (the percentage of total variability due to heterogeneity; Higgins and Thompson, 2002). All calculations were conducted in R (Version 2.15.1, R Core Team, 2012). A sufficient amount of data of the combined effect of OA and temperature within global literature review allowed for a similar meta-analysis at co-occurring temperature and Ω_ar levels, but not for DO. To account for an additional warming effect, we have conducted the same steps as described above for OA meta-analyses (Table 2). In short, we first calculated OA effect independently under ambient temperature, and then compared it with the OA effect under an elevated temperature.

Table 1. Responses considered for threshold analysis.

Table 2. Temperatures and response measures of studies used for meta-analysis. Only L. helicina studies have been included.

Breakpoint Analyses

Pteropods experience variable exposure to low Ω_ar and low DO and occasional thermal stress in their natural environment. As used here, a threshold is defined as the point at which there is a statistically significant change in specific pteropod response (Y-axis) for an incremental change in the environmental stressor (x-axis; e.g., Ω_ar). The goal of breakpoint analyses was to identify thresholds from each set of published studies, documenting Ω_ar stress effects on pteropod responses. To link pteropod responses with Ω_ar stress specifically, we have excluded the remaining responses that showed either insignificant response in meta-analyses, with datapoints from these studies not being included in the subsequent breakpoint analyses. We only focused on the Ω_ar stress-related threshold derived from negatively affected responses, but throughout the study, we acknowledged that two insignificant results also exist (e.g., excretion, respiration).

For experimental data, we evaluated the transition from no observed effect (experimental control) to statistically significant least observed effects (LOE) and severe observed effect (SOE, Supplementary Figure 1). For field-based stress-response studies, piecewise regression analysis was used to identify the breakpoint in slopes in the response measures over the gradient in OA stress (package “segmented”; version 2.15.1, Muggeo, 2008; R Core Team, 2012). For the same measures, we retained information on duration during threshold derivation for each pathway all thresholds were duration dependent. While meta-analysis was used to study the impact of confounding effects on the threshold, the result of meta-analysis was not directly applied to the breakpoint analyses.

Expert Assessment of Uncertainty and Determination of Appropriate Exposure Duration

An expert working group, with their contributions acknowledged through their co-authorship in this paper, was assembled at the outset of the synthesis process. The scientists that comprised this group were chosen based on depth of knowledge and expertise in pteropod physiology and responses to OA and multiple stressors, as well as their demography, population dynamics, and life history strategies. The expert working group added unpublished studies to the database and provided a review process for both the meta-analysis and threshold analyses (e.g., Teck et al., 2010). This type of check was important because all of the studies were treated equally in the analyses, despite potentially unequal levels of data quality and quantity. The experts were provided with the opportunity to refine the interpretation of the analytical outcomes based on their knowledge of how those studies should be weighted. For each response measure, each expert identified the OA threshold they would assign based on the breakpoint analysis. In cases where the experts disagreed with the analytical result and were divided (which was rare) the threshold value with the most votes was adopted; if there was not a clear majority, the value was discussed amongst the experts until a consensus was reached. Meta-analysis greatly simplified the work of experts, enabling the identification of thresholds from 22 response measures and duration combinations down to five major endpoints represented on the final graphic.

In addition, there was an unequal quantity of study results, as well as differing levels of agreement among studies, for each of the pathways examined. To enhance the applicability for each of the findings, the experts were asked to provide a confidence level for each pathway threshold. This was done using the Intergovernmental Panel on Climate Change (IPCC) confidence model. This metric qualitatively expresses the confidence in the validity of a finding, on the type, amount, quality, and consistency of evidence to assign a confidence level (Supplementary Figure 2). An uncertainty rating was assigned based on IPCC uncertainty guidance that evaluates the confidence of findings in terms of evidence and agreement (low, medium, high; Mastrandrea et al., 2010). Meta-analysis results were used as an aid to discuss the confidence score.

Finally, the experts were asked to identify the most appropriate time frame over which the thresholds should be applied. The range of exposure durations varied among the underlying studies and expert judgment was needed to determine whether these were inherently acute as opposed to chronic exposure responses, and the extent to which fluctuations above and below these thresholds were cumulative.

Applying Thresholds to Biogeochemical Model Outputs

Experts discussed guidance that should be provided in order to achieve a more consistent application of these threshold to ocean observations and ocean numerical model output. They recommended use of measures such as intensity, duration, and severity of departure from the thresholds, based on Hauri et al. (2013). Intensity represents the magnitude of the departure from the threshold, calculated as the negative average deviation from the threshold. The duration represents the number of days below the threshold expressed as a % of the total analysis period. The severity is the product between the intensity and duration, with the % time of Ω_ar below the threshold deviation. Calculations of duration, intensity, and severity of departure from the derived OA thresholds were carried out through the use of a set of model simulations from the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005, 2009), with biogeochemical elements (Biogeochemical Elemental Cycling; BEC, Moore et al., 2004) that has been developed for the CCS, the Pacific Coast of North America. The thresholds were particularly applied to the southern part of the CCS, particularly for the region of California.

The ROMS-BEC model provides a realistic three-dimensional representation of the physical circulation (Renault et al., n.d.) and the biogeochemical cycles of nutrients, oxygen, inorganic, and organic carbon species, and several plankton functional groups (Deutsch et al., n.d.). Ω_ar was calculated based on temperature and oxygen fields of the ROMS-BEC model using the Juranek et al. (2009) algorithm. In separate publications, we describe the model set up, and demonstrate that it captures the broad patterns of physical circulation, productivity, dissolved nutrients, and DO along the CCS, and their temporal variability and trends (Deutsch et al., n.d.; Renault et al., n.d.). The simulation analyzed here is a downscaling of a 4 km-resolution US West Coast solution to the southern portion of the California coast (Kessouri et al., n.d.), with data deriving directly from the models. It runs with a spatial resolution of 1 km, submesoscale permitting, and spans the 10-year period from 1997 to 2006 (Deutsch et al., n.d.; Renault et al., n.d.) at 1 km spatial resolution and 1-day average outputs (Kessouri et al., n.d.). Daily Ω_ar averages of 0–200 m depth in the water column were used either over 2 periods for mild dissolution and survival, respectively over the 10-year period from 1996 to 2007. The threshold diagnostics included two case studies. First, juvenile mild dissolution (Ω_ar threshold of 1.5 for a duration of 5 days) was applied for a 2-month period, between April and May, which coincides with the presence of early life stages (larval and juvenile) along the southern part of CCS (Wang et al., 2017) that are most susceptible to mild dissolution, for a total of 610 days over the 10-year simulation. Second, threshold for severe dissolution (Ω_ar threshold of 1.2 for 14 days duration) was applied to the averaged 0–200 m water column over the same time period and across the same spatial scales as the mild dissolution threshold.

Results

Omega Saturation State as the Consensus Measure of OA Stress

All 15 studies selected for the final threshold identification used aragonite saturation state (Ω_ar) as the chemical measure of organismal stress (Figure 2). The majority of focal endpoints had demonstrated response to Ω_ar, which was also a preferred proxy because the pH and pCO₂ only showed unambiguous responses for a few sublethal endpoints (e.g., metabolic and behavioral responses). Ω_ar response in experimental treatments spanned a range of under- and super-saturated conditions (Ω_ar of 0.7–3.2; Figure 2A), with exposure duration varying from days to months (Figure 2B).

Meta-Analysis Informs Evaluation of Response Measures and Evaluation of Confounding Factors

After pooling pteropod data across species, life stages, and geographic locations, a significant negative response to Ω_ar was found across most exposure responses, physiological and lethal processes, indicating a high signal-to-noise ratio (Figure 3). Severe shell dissolution, egg developmental impairment, larval growth, and survival had the greatest magnitude of Logarithmic Response Ratio (LnRR) response to Ω_ar, and thus were prioritized for threshold development. In contrast, a few processes had insignificant LnRR, interpreted as either poor signal-to-noise ratio (e.g., excretion rates), or insufficient datapoints (e.g., larval survival), and were excluded from threshold synthesis. Finally, some processes showed either a significant positive response of the stress (gut clearance), or insignificant LnRR response (respiration, excretion rate). The threshold synthesis focused on deriving thresholds only for significant negative biological responses. Our results for pteropods show great uniformity of within-taxa impairment responses to OA, lending support for the decision to pool data across experiments when considering the thresholds during the final synthesis phase. We elaborate on these four aspects of response below.

Figure 3. (A) Forest plot of Ω_ar effect on various biological processes across exposure response, sublethal, and lethal responses in pteropods. The numbers in the parenthesis indicate the number of experiments for each response measure, error bars indicate 95% confidence intervals. LnRR (Logarithmic response ratio) indicates the magnitude of response with its significance (p < 0.05) indicated by the asterisks. Percentage values along the top of the figure indicate percent decrease in the given endpoint, compared to the control. (B) “Effect” (box-and-whiskers) and “control” (rectangles) conditions as used in the experimental conditions, with the dotted line positioned at Ω_ar = 1 and 1.8 for comparison with known biogeochemical Ω_ar thresholds.

The magnitude of sublethal and lethal responses was predominantly life-stage specific. Larvae and juveniles were substantially more affected by Ω_ar than adults, while egg sensitivity was comparable to adults (Figure 4). We observed the greatest reduction in survival in juveniles, and significant for larvae, whereas no effects were detected for adults, consistent with Dupont et al. (2010) and Kroeker et al. (2013). A pronounced sensitivity of larvae and juveniles suggest that the Ω_ar exposure of early pteropod life stages could represent a population bottleneck, a key consideration to translate effects from individuals to population level (sensu Bednaršek et al., 2016). For this reason, threshold analyses were life stage-specific for individual endpoints, when meta-analysis could not otherwise justify their combination.

Figure 4. The variability in LnRR responses across (A), exposure response, sublethal, and lethal responses, (B), across various life-stages, and (C), with specific regionality. LnRR (Logarithmic response ratios) indicate the magnitude of the response, with the asterisks indicating its significance (negative associated with negative LnRR values; and positive associated with positive LnRR values). LnRR value crossing zero line indicates the insignificance of the response.

Species-specific responses represent additional sources of variability, as demonstrated by highly variable LnRR, albeit with an overall significant negative response, except for significantly less sensitive temperate species (L. retroversa; Figure 4). Species found in the upwelling and polar regions had the greatest overall OA response, indicating similar vulnerability across habitats with some of the most severe OA gradients. In addition, tropical species also show a significant negative sublethal response, comparable with an LnRR to their west coast counterparts (Figure 4). The LnRR of temperate species was only significant when data on highly variable respiration responses were removed from the meta-analysis. LnRR differences on lethal endpoints between the species of divergent geographic origin suggest local populations may have a differential response to Ω_ar stress. Despite this, significantly negative LnRR responses across species provided evidence to combine experimental studies during the threshold synthesis stage.

Finally, meta-analysis of six different response measures combining thermal with Ω_ar stress significantly increased the magnitude of severe dissolution and mortality (Figure 5), changed the direction of the LnRR for the sublethal responses (e.g., increased respiration), or resulted in no significant effect (e.g., calcification, excretion growth). Meta-analyses cannot distinguish between these additive and synergistic effects, but evidence exists that these two stressors can synergistically impact sublethal and lethal endpoints (Lischka and Riebesell, 2012).

Figure 5. Meta-analysis of OA (Ω_ar) and warming effects on the variety of response measures at ambient and elevated temperature, with underlying Ω_ar conditions (A) Forest plot of Ω_ar and thermal effect across exposure response, sublethal, and lethal responses in pteropods. The numbers in the parenthesis indicate the number of experiments for each response measure. LnRR (Logarithmic response ratios) indicates the magnitude of response with its significance (p < 0.05) indicated by the asterisks. (B) “Effect” (box-and-whiskers) and “control” (rectangles) conditions as used in the experimental conditions, with the dotted line positioned at Ω_ar = 1 and 1.8 for comparison with known biogeochemical Ω_ar thresholds.

Synthesis of Thresholds by Endpoints

Expert synthesis identified duration-dependent thresholds for Ω_ar effects on dissolution, calcification, growth, egg development, and survival (Figure 6), ranging from 1.5–0.9, with assigned uncertainties (Figure 7; Supplementary Figure 2).

Figure 6. Synthesis of thresholds related to six different endpoints, characterized by its magnitude (Ω_ar), and duration of exposure response. The thresholds were determined based on experimental and field data using breakpoint analyses and expert consensus.

Figure 7. Confidence score based on the combination of evidence and agreement for six different endpoints as determined by the expert consensus. Figure adapted from IPCC 2014 report.

Mild and Severe Shell Dissolution

Five different exposure measures were combined into either mild (including opacity, surface shell dissolution and transmittance) or severe (deeper) shell dissolution (Table 1). Based on comparable magnitude of LnRR across disparate measures of dissolution, life stages and species, meta-analysis provided support for aggregating various measures. The presence of mild shell dissolution is the most sensitive measure, providing an early warning response, while severe dissolution imposes sublethal impairment because it likely requires additional energetic costs for repair and maintenance (Lischka et al., 2011; Bednaršek et al., 2016; Peck et al., 2018).

Mild and severe shell dissolution processes differ substantially in their thresholds, which ranged from 1.06 < Ω_ar < 1.59 for 2–8 weeks, respectively (Bednaršek et al., 2014a; Bergan et al., 2017; Table 1). These experimentally derived thresholds were independently verified from a field stress-response study, with derived changepoint threshold of Ω_ar = 1.3 (Bednaršek et al., 2014b; Supplementary Figure 1). There was robust consensus among the experts with high levels of agreement, that dissolution is not species- or life stage-dependent (Figure 7). The thresholds for mild and severe dissolution were set at Ω_ar = 1.50 for 5 days and Ω_ar = 1.20 for 14 days, respectively, with the highest confidence score (99%) assigned (Figure 7).

Calcification

In comparison with dissolution, calcification response was less sensitive and less consistent to Ω_ar, the latter attributable to differences in life history and energetic state characteristics. This implies variable control over calcification processes, likely because calcification is energetically expensive when sustained at low Ω_ar gradients (Lischka et al., 2011; Bednaršek et al., 2017). While juveniles showed reduced calcification at low Ω_ar (Comeau et al., 2010a; Bednaršek et al., 2017), an increase in calcification was seen in the less sensitive adults of L. retroversa (Maas et al., 2018). Experimentally derived LOE ranged from Ω_ar = 1.25 for 3 days (Comeau et al., 2010a) to Ω_ar = 1.02 for 1-3-week duration (Maas et al., 2018), which coincides with a significant calcification decline also in Heliconoides inflatus (Moya et al., 2016). These duration-dependent thresholds were validated with a field-derived changepoint of Ω_ar = 1.02 (SE = 0.07; Bednaršek et al., 2017). While calcification processes might still be confounded by the species- and life stage-responses, meta-analysis supported a decision that the thresholds could be used across potentially confounding parameters. As such, calcification was positioned at Ω_ar = 1.25 for 7 days of exposure. This threshold varied from medium-to-robust evidence and medium-to-high agreement, with final expert consensus giving it an overall medium-high confidence score.

Growth Reduction

Shell growth parameters in juveniles included shell diameter, length, or their ratios (Lischka et al., 2011; Lischka and Riebesell, 2012; respectively), while the larval growth was measured using calcein (Comeau et al., 2010b). Shell growth was negatively correlated to Ω_ar and showed a uniform response; meta-analysis demonstrated that responses were confounded by life stage but not by species. Larval growth impairment showed a significant difference in the magnitude of response in comparison with adults and was one of the most sensitive to Ω_ar. Of those studies measuring Ω_ar growth impacts, data could only be obtained for the study by Manno et al. (2012) but not from Thabet et al. (2015) and Comeau et al. (2010b) because LOE could not be extracted, or the uncertainty was too great, yet both studies demonstrated a significant decrease in larval growth and developmental delay. The derived LOE based on Manno et al. (2012) was Ω_ar = 1.15 for 7 days. In addition, field and laboratory studies agree with reduced shell length and thickness at Ω_ar ~ 1 (Lischka et al., 2011; Lischka and Riebesell, 2012; Bednaršek et al., 2017), providing a solid cross-validation of the expert-selected threshold. The growth threshold is readily applicable for slower-growing polar species but might be more difficult to extrapolate at different latitudes, given species-specific growth patterns. An additional uncertainty is attributable to disparate sizes of experimental organisms among the studies. Therefore, the growth threshold was assigned medium evidence and medium-to-high agreement, yielding an overall medium-to-high confidence score. While the uncertainty was too large to identify a threshold specifically for larval growth, we recognized its importance as an independent endpoint.

Severe dissolution, calcification, and growth showed close proximity of the Ω_ar thresholds, all occurring within ±0.1 unit and within a duration of 7–14 days. This proximity of the thresholds indicated that all of these processes are either simultaneously affected or one process results in a cascade of responses, provided sufficient duration of exposure. For example, while calcification is a more sensitive process, it can affect growth through a trade-off, such that the organisms' use of additional energy for calcification results in a reduced energy scope for growth. The threshold for mild dissolution is much higher than for severe dissolution, indicating this to be predominantly a geochemical process that is initiated by changes in Ω_ar in external conditions rather than organismal response.

Pre-larval Development

Maternal OA exposure during only the egg-brood phase resulted in smaller eggs with lower carbon content, while embryonic OA exposure during the embryonic development phase retarded the development rate. Combined maternal and embryonic exposure to Ω_ar < 0.90 for 3–5 days reduced the percentage of eggs successfully reaching the organogenesis stage by 80% (Manno et al., 2016). Although this study demonstrated the effect of maternal exposure on egg development, the threshold was at Ω_ar < 0.90 for 2 days. Given the limited evidence, the overall confidence score for this threshold was medium-to-low. As such, use of the egg development threshold is caveated, particularly with respect to species difference; it is relevant for species (e.g., L. helicina) in which shell development occurs within the egg capsule, but not for other species (e.g., L. retroversa) that hatch to unshelled larvae after spawning and develop into veligers with a cup-shaped shell or are egg brooders with a protected egg development phase, as in the case of H. inflatus (Lalli and Gilmer, 1989).

Survival

The final consideration of survival thresholds was based on L. helicina juveniles from polar and upwelling regions (Lischka et al., 2011; Bednaršek et al., 2017). The LOEs ranged from Ω_ar = 0.94 for a 4-week lab experiment (Lischka et al., 2011) to Ω_ar = 1.06 for 6-week exposure in a field experiment (Bednaršek et al., 2017). The juvenile and adult survival threshold was set to Ω_ar = 0.95 for 14 days. Meta-analysis results showed significant variation of LnRR in survival between the juveniles and adults, as well as for the species-specific response with greater lethality in upwelling species.

Although several studies examined the survival of pteropods, the range of responses varied. No effect on survival was found when the experimental duration was <1 week (Seibel, unpublished data; Manno et al., 2012). Survival rates across studies depended on the difference in nutritional and energetic status of organisms (Seibel et al., 2012) and their history of exposure (Bednaršek et al., 2017). This affected the intercept on the Y-axis (% survival) but neither the slope of the relationship nor the threshold, which was comparable for both studies. In addition, a significant decrease in larval survival was reported in only one study (Thabet et al., 2015), but due to the high variability of this study, it was excluded from the threshold synthesis. Because survival responses are life-stage specific and attributed to a variety of confounding parameters, such as food availability, temperature, or exposure history, and a well-documented captivity effect of cultured pteropods, the survival threshold was designated a low confidence score.

Demonstration of Threshold Application to Ocean Model Numerical Simulations

Application of the duration, intensity, and severity of departure of the mild dissolution and survival thresholds to ROMS-BEC numerical simulations of the CCS illustrates key points in their application (Figure 8). First, a comparison of thresholds of mild (Figures 8A–C) vs. severe dissolution (Figures 8D–F) applied to modeled Ω_ar represents bookends of the range of impacts related to shell processes. While conditions appear to exist that frequently induces mild dissolution in pteropod vertical habitat, the incidence of crossing severe dissolution threshold was less common, mostly along the onshore region at a lower part of the 200 m vertical habitat. Second, each panel provides spatially explicit information on the magnitude of departure from the threshold (Figures 8A,D), the percent of time it occurs (Figures 8A,E) and the severity of those departures, combining the magnitude and duration of exposure into one metric (Figures 8C,F).

Figure 8. Examples of threshold applications to modeled estimates of Ω_ar at 200 m off the entire California coast. Magnitude, duration, and severity of exposure were calculated for two thresholds, respectively: (A–C), Mild dissolution threshold of juvenile pteropods (Ω_ar = 1.5; applied for 5 days over the period of March to May); (D–F), severe dissolution of juvenile pteropods (Ω_ar = 1.2; applied for 14 days over the period of March to May).

As illustrated here, there is a consistent pattern of increasing magnitude, duration, and severity from offshore to onshore. For both processes (mild and severe dissolution), the threshold represented on the 2-dimensional maps of the saturation state of Ω_ar show deviations from the threshold in the magnitude, duration, and severity of the seasonal exposure that coincided with pteropod spawning. The model shows that, on average, the greatest magnitude, and severity of exposure is in the shelf regions, with approximately 80 and 20% time exposure below the threshold; and 30 and 20% deviation below the threshold for the mild and severe dissolution, respectively. This implies that the early life stages are most susceptible to dissolution during their most vulnerable periods.

Discussion

Global Synthesis of the Effects of OA on Pteropods: Synthesis and Comparison With Other Organism Groups

Meta-analysis and global literature review of threshold demonstrated the high sensitivity of pteropods across multiple responses to OA. The meta-analysis demonstrated (1) significant, negative responses to decreasing Ω_ar across most response measures, (2) the magnitude of response was duration and life-stage specific, and (3) temperature enhanced the magnitude of two responses negatively impacted by Ω_ar, shell dissolution, and survival. Comparison of pteropod meta-analyses across broader marine organisms (Kroeker et al., 2013) shows that pteropods bore very similar responses with most calcifying groups, including negative responses across survival, calcification, growth, and development. The greatest comparable magnitude of responses were with the diverse group across the entire mollusk group, especially the sensitivity of the early-life stages. In comparison with the other groups, pteropods show more sensitivity in developmental response. When temperature responses are included in the meta-analysis, the trend is noticeable toward lower survival in both pteropods (this study) and the other calcifiers (Kroeker et al., 2013); however, the difference between OA effects of control and elevated temperature does not explain the significant amount of variation in those responses.

With respect to the threshold comparison with previous studies, oysters (Pacific, Olympia and Eastern oyster) and mussels showed great similarity in the magnitude of thresholds with our current study, although the exposure duration for these species have not been defined. Pteropod mild shell dissolution at Ω_ar = 1.5, an early warning response, is comparable to sublethal thresholds for oysters or mussel larvae in the Ω_ar range of 1.4 and 2 (Barton et al., 2012; Hettinger et al., 2012) and Ω_ar = 1.8 (Gaylord et al., 2011). In terms of acute responses, the pteropod survival threshold of Ω_ar = 0.95 was lower (i.e., occurring at higher Ω_ar stress), yet analogous to bivalve larvae sublethal acute threshold at Ω_ar = 1.2 (Waldbusser et al., 2015). The differences are life-stage specific; larval oyster stages have the greatest sensitivity, while pteropod thresholds reflects a sensitivity across all life stages. Higher threshold values (Ω_ar > 1.5) are referenced in consideration with the cumulative effects of all biological processes that can compromise the overall biological status (Gimenez et al., 2018). As such, our pteropod thresholds are valid, yet conservative estimates for individual endpoints, but do not address the cumulative impact of impairment pathways.

Advantages to Applying Thresholds to Biogeochemical Model Output

Scientific consensus on OA thresholds is key for their application to OA status and trend assessments, predictions of the effects of climate change, assessment of impact of local pollution inputs, proposed remediation strategies, and as the foundation for OA water quality goals (Weisberg et al., 2016). However, guidance is required for their consistent application, particularly with respect to decisions on how data should be aggregated to assess the spatial extent and temporal duration needed to trigger a biological endpoint. Numerical model simulations were used to illustrate the basic metrics of duration-dependent threshold applications (magnitude, duration, severity; Hauri et al., 2013) and can be used to indicate reduced habitat suitability for pteropods. Appropriate vertical aggregation of observations or modeled output remains a key question. In our particular case, experts delineated the maximum depth limit to which these thresholds should be applied (<500 m), noting that the most abundant pteropod species are typically limited to the upper 200 m of the water column (Bednaršek et al., 2012a); however, the experts were challenged to provide simple guidance on what statistical approaches to use when aggregating the data across the vertical or spatial scales to demonstrate the reduction in suitable habitat that would result in population level effects. The challenge lies in a variety of complex factors that can impact population response that work across time scales not necessarily aligned with the observed responses.

One of the most important considerations is to understand these thresholds within the context of the pteropod life cycle, particularly the timing of reproduction and growth against the seasonal Ω_ar settings, related to upwelling cycles (Feely et al., 2008), primary productivity, freshwater runoff patterns, and ice melt (Harris et al., 2013). Thresholds for survival, growth, and egg development are life stage specific, while dissolution and calcification can be applied across life stages. Different life stages can withstand different levels of mortality without significant long-term impacts to the population (Batiuk et al., 2009), so identification of the critical conditions under which population-level effects are predicted will depend on stage sensitivities and cumulative effects of co-occurring environmental conditions. For example, the reproduction of L. helicina in the northeast Pacific occurs in the spring and fall (Wang et al., 2017). With the early-life stages having the highest sensitivity to survival and development, exceeding the Ω_ar threshold during these periods will determine the population outcomes. In addition, a prolonged duration of exposure across several life stages might first carry-over through maternal effects toward higher juvenile vulnerability, followed by cumulative effects throughout the life history to lead toward the population effects, particularly where pteropod reproduction is continuous throughout much of the year, e.g., in the Southern Ocean (Bednaršek et al., 2012b). Demographics modeling can help identify the most critical life stages in a population (e.g., Bednaršek et al., 2016) that, in combination with the predicted seasonal changes in Ω_ar and their regional variability (Kwiatkowski and Orr, 2018), help identify critical timing and regional aspects of the population bottlenecks. It can also be used to develop more precise recommendations for the vertical aggregation of output or observation or spatial extent needed to estimate a population impact.

Limits of This Synthesis and of Threshold Analysis

On several occasions, experts did not support assigning the thresholds to a specific process, or thresholds had high uncertainty. These instances were due to insufficient or poorly-curated data, disparity in response measures or equivocal responses with a poor signal-to-noise ratio, or data gaps that did not allow for the extrapolation across spatial scales or life stages. Furthermore, interpreting the large variability in responses (poor signal-to-noise ratio) proved difficult. An example of this is the variation in metabolic rate; several confounding parameters could have impacted metabolism processes, yet they were not sufficiently defined or understood during the experimental stage, such as the nutritional or energetic status. As such, consideration of respiration thresholds was eliminated from the final synthesis (Comeau et al., 2009; Moya et al., 2016; Lischka and Riebesell, 2017). We also recognized some additional caveats from the experimental studies that need to be considered; for example, experimental datasets were skewed toward shorter experiments, which made inferences of chronic effects less certain, or insufficient feeding that might have generated more acute responses that would have otherwise occurred in their natural environment. Finally, an insufficient number of experiments have explicitly tested multiple stressor effects on pteropods (e.g., Lischka and Riebesell, 2012, 2017; Manno et al., 2012; Maas et al., 2018), making the stated thresholds only representative of Ω_ar in a single stressor future ocean. Approaches such as regional-specific habitat suitability indices that consider multiple stressors may represent a useful complementary approach. Experts agreed that currently synthesized thresholds are globally applicable to calcifying pteropods across varying spatial scales. Additional field and experimental data can refine these thresholds over time, particularly if conducted in a way to extract or identify thresholds rather than simply document OA adverse effects.

Currently, insufficient experimental data exist to characterize DO and thermal stress effects on OA thresholds. Several thermal and OA stress studies allowed for meta-analysis to be conducted as a first step but that was not feasible for an OA- low DO combination. Furthermore, the breakpoint analyses were not feasible for Ω_ar in combination with either thermal or low oxygen data. The stress response is expected to occur at lower stress under multiple stress conditions that will result in higher threshold values. It is extremely important to account for all possible co-occurring stressors in the species' natural habitat to yield accurate predictions of impaired biological processes. The single Ω_ar thresholds generated in this study are likely biased with respect to the less detrimental responses and will need to be upgraded within the multiple stressor context when more data is available. However, a few field studies of pteropod vulnerability are available to provide insights into potential changes of the threshold sensitivity due to multiple stressor effects. Under combined low DO and OA, some species showed metabolic depression (Maas et al., 2016), while this was a prevalent response for all the species under low DO and warming (Maas et al., 2012). In addition, increased pteropod mortality and shell dissolution was demonstrated for pteropods under the combined effect of warming under marine heat wave and OA. OA stress in combination with low salinity also impacted pteropod behavior (Manno et al., 2012). Based on these studies, potentially higher threshold values indicative of induced stress at less unfavorable conditions under combined stressors are expected for the most sensitive processes, such as dissolution and survival, and possibly neurologically induced behavioral impairments that is already negatively affected under OA alone (cf. Moya et al., 2016; Bergan et al., 2017).

These kinds of integrated analyses provide an opportunity for identifying information gaps; here we identified three areas where future research investments would substantially improve our ability to define OA thresholds. First, most of the experiments measure effects based on static Ω_ar conditions, which are not representative of realistic conditions. As such, variable exposure studies are needed to evaluate the impact of Ω_ar natural variability, especially because there is a complete lack of understanding of how dynamic Ω_ar exposure effects the level of resiliency. Resiliency may occur if the relaxation periods between stressor exposure events are of sufficient length to allow partial recuperation before the next exposure event, thereby lessening its susceptibility. Second, since the bulk of the studies were experiments conducted on individual organisms, thresholds derived from field-based stress-response observations are needed to validate the relationship between OA exposure and biological endpoint or population measures (e.g., percent individuals affected, abundance, etc.). This should be tested across life stages and species found at various latitudes to validate the applicability of thresholds on a global scale. Third, process studies and population modeling are needed to provide an integrated view of how the physiological responses to environmental stressors, coupled with seasonality exposure, translate to population-level effects. These studies are central to developing a solid scientific foundation for multi-trophic level ecosystem models to support investigations of climate change related ecological tipping points.

Author Contributions

NB conducted the literature review, statistical analyses, and data syntheses. NB and NN conducted meta-analyses. FK performed modeling and provided model outputs. NB led manuscript writing, with the other co-authors RF, EH, BH, PL, SL, AM, KM, MS, and SW contributing to writing and editing.

Funding

This research was supported by California Ocean Protection Council, grant number C0302500. RF was supported by the NOAA Ocean Acidification Program.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dana Shultz for her assistance in manuscript preparation. NOAA PMEL Contribution number is 4969.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2019.00227/full#supplementary-material

References