Integrated Assessment of Ocean Acidification Risks to Pteropods in the Northern High Latitudes: Regional Comparison of Exposure, Sensitivity and Adaptive Capacity

Abstract

Exposure to the impact of ocean acidification (OA) is increasing in high-latitudinal productive habitats. Pelagic calcifying snails (pteropods), a significant component of the diet of economically important fish, are found in high abundance in these regions. Pteropods have thin shells that readily dissolve at low aragonite saturation state (Ω_ar), making them susceptible to OA. Here, we conducted a first integrated risk assessment for pteropods in the Eastern Pacific subpolar gyre, the Gulf of Alaska (GoA), Bering Sea, and Amundsen Gulf. We determined the risk for pteropod populations by integrating measures of OA exposure, biological sensitivity, and resilience. Exposure was based on physical-chemical hydrographic observations and regional biogeochemical model outputs, delineating seasonal and decadal changes in carbonate chemistry conditions. Biological sensitivity was based on pteropod morphometrics and shell-building processes, including shell dissolution, density and thickness. Resilience and adaptive capacity were based on species diversity and spatial connectivity, derived from the particle tracking modeling. Extensive shell dissolution was found in the central and western part of the subpolar gyre, parts of the Bering Sea, and Amundsen Gulf. We identified two distinct morphotypes: L. helicina helicina and L. helicina pacifica, with high-spired and flatter shells, respectively. Despite the presence of different morphotypes, genetic analyses based on mitochondrial haplotypes identified a single species, without differentiation between the morphological forms, coinciding with evidence of widespread spatial connectivity. We found that shell morphometric characteristics depends on omega saturation state (Ω_ar); under Ω_ar decline, pteropods build flatter and thicker shells, which is indicative of a certain level of phenotypic plasticity. An integrated risk evaluation based on multiple approaches assumes a high risk for pteropod population persistence with intensification of OA in the high latitude eastern North Pacific because of their known vulnerability, along with limited evidence of species diversity despite their connectivity and our current lack of sufficient knowledge of their adaptive capacity. Such a comprehensive understanding would permit improved prediction of ecosystem change relevant to effective fisheries resource management, as well as a more robust foundation for monitoring ecosystem health and investigating OA impacts in high-latitudinal habitats.

Introduction

The pervasiveness of global climate change threatens biodiversity, habitat suitability, ecosystem function, structure and services (Bindoff et al., 2019). Such changes are especially pertinent in high-latitudinal habitats across the eastern North Pacific Ocean, including the Gulf of Alaska (GoA) and its subpolar gyre, the polar waters of the Bering Sea, one of the most productive seas on the planet (Brown et al., 2011), and the Amundsen Gulf of the Canadian Beaufort Sea. Yet, the subpolar and polar seas are sensitive to ocean acidification (OA) due to naturally low buffering capacity, coupled with multiple climate-change and related physical-chemical drivers, such as increasing freshwater input from glacial runoff, rivers, sea ice melt, and increasing temperature (Mathis et al., 2011; Evans et al., 2013; Cross et al., 2013, Hauri et al., 2020). Substantial changes in carbonate chemistry conditions have already been detected, demonstrating increasing corrosiveness in high-latitudinal North Pacific regions (IPCC, 2014, 2019; Carter et al., 2017, 2019; Bindoff et al., 2019; Hauri et al., 2021). These changes coincide spatially and temporally with the habitats of important commercial and subsistence fisheries (Mathis et al., 2014), underscoring the need for a detailed understanding of the species tolerances and their habitats that will most likely be at risk. Currently, experimental data from the Pacific subpolar gyre, GoA, and Bering Sea show that economically important crab species (Long et al., 2013, 2016) and coho salmon (Williams et al., 2019) are sensitive to OA. However, very limited information exists from field measurements on the current impacts of OA on the most sensitive marine calcifying zooplankton (Cai et al., 2020), such as pteropods, with direct or indirect importance for the fisheries. Pteropods are a key species within the high-latitudinal ecosystems (Kobayashi, 1974; Fabry, 1989; Bednaršek et al., 2012a) that demonstrate severe sensitivity to OA. Multiple biological pathways of this species are impacted by OA conditions, compromising their growth, fecundity, and ultimately their survival (Comeau et al., 2009; 2010; Lischka et al., 2011; Bednaršek et al., 2014, 2019; Feely et al., 2016; Manno et al., 2016). As a result of such OA-related sensitivity, pteropods rank in the top quartile of ecological integrity indicators and are used to identify early warning signs and habitats of concern related to OA (Bednaršek et al., 2017).

Pteropods have a wide-ranging, year-around distribution in the subpolar gyre (Mackas and Galbraith, 2012), offshore and coastal habitats of the GoA (Tsurumi et al., 2005; Doubleday and Hopcroft, 2015), and Bering Sea (Eisner et al., 2014, 2017; Janssen et al., 2019), and the Amundsen Gulf (Darnis et al., 2008; Walkusz et al., 2013; Niemi et al., 2021). Pteropods in the GoA are commonly found in the diets of various life stages of Pacific Salmon species, such as pink, sockeye and chum salmon (Auburn and Ignell, 2000; Daly et al., 2019), that rely upon pteropods to support somatic growth and lipid reserves, especially during the critical summer and autumn growth periods (Auburn and Ignell, 2000; Armstrong et al., 2005; Beauchamp et al., 2007; Zavolokin et al., 2014; Doubleday and Hopcroft, 2015; Daly et al., 2019). They can also play a role in the diet of various fish in the Bering Sea, where high pteropod abundance occurs in the southern part of the basin, predominantly during the autumn (Eisner et al., 2014, 2017, 2017), making them one of the most dominant components of the diet of various salmon species (Zavolokin et al., 2007). In the Amundsen Gulf, pteropods comprise a small component of demersal Arctic cod diets (Majewski et al., 2016) and can be found at high numbers in the stomachs of coastal fish, including Arctic char (Niemi et al., 2021). In addition, given that pteropods in these regions have one of the highest abundances and biomass on the global scale, they also significantly contribute to the biological pump, especially to the carbonate flux (Bednaršek et al., 2012a). The most abundant species is Limacina helicina, which is present in two forms that are morphologically distinct from each other in the adult stage; one designated as L. helicina helicina, distinguished by an extended spire with striated surface, and the other as L. helicina pacifica, with a depressed (flat) spire and no striae present (McGowan, 1963; Janssen et al., 2019).

Assessing OA risks requires explicit linkages between three components (sensuWilliams et al., 2008): (1) regional and local exposure to OA (e.g., changes in habitat suitability, spatial and temporal habitat change); (2) species’ sensitivity governed by intrinsic metabolic and ecological traits, such as climatic preferences, reproductive outputs, habitat use, and biotic and abiotic interactions; and (3) resilience and adaptive capacity (e.g., genetic and phylogenetic diversity, phenotypic plasticity, life history traits, dispersal potential, spatial connectivity). Emerging literature illustrates the species’ sensitivity under increasing OA exposure, but little is known about a vital component of the risk assessment: their distinct ecological, phenotypic or evolutionary specialization to changing conditions (e.g., Williams et al., 2008; Burridge et al., 2015).

Given the key role of pteropods in high-latitude ecosystems, there is an urgent need for an integrated assessment of OA risks for pteropod population viability in the high-latitudinal habitats, including the eastern part of the North Pacific, GoA, Bearing and Amundsen Gulf. Following Williams et al. (2008), we delineated pteropod OA risk assessment as a combination of exposure, sensitivity, and adaptive capacity. First, OA exposure was estimated using hydrographic observations and regional ocean biogeochemical models; we used model outputs to define short-term, seasonal OA exposure, as well as longer-term OA trends with the interannual and decadal changes that characterize a gradual OA intensification in the species’ exposure. These metrics of change offer tools to delineate the location of OA hotspots and habitat suitability with respect to OA, as well as temporal windows during which the biological exposure will be most significant.

Second, sensitivity was determined by observing early and well-defined negative effects of OA on pteropods, starting at aragonite saturation state (Ω_ar) values below 1.3 (Bednaršek et al., 2019), which is related to the shell biomineralization processes (dissolution and calcification). Such sensitivity was explored with respect to various morphotypes and matched with Ω_ar observations of the pteropods’ place of origin. Another important aspect of sensitivity is life history, specifically that of exposure severity as it co-occurs with varying levels of vulnerability over pteropod life stages.

Third, potential resilience was related to adaptive capacity, or phenotypic plasticity, genetic population structure, and population (spatial) connectivity. Genetic-based analyses can be used to identify cryptic species and their distribution and connectivity across all regions. Given limited knowledge about pteropod cryptic diversity (on the level of Limacina spp.) in the region, we focus our initial efforts on characterizing diversity at the species, subspecies, and morphotype levels (Maas et al., 2013; Gasca and Janssen, 2014; Sromek et al., 2015a,b; Shimizu et al., 2018, 2021). OA-induced exposure would not necessarily result in local vulnerability if individuals were replaceable. Evidence of co-occurring cryptic species or subspecies with differential response to OA might buffer changes in abundance. Similarly, large-scale spatial connectivity between geographic areas may promote replacement from unaffected areas.

Here, we conduct a comprehensive assessment of OA risks for the pteropod population in the eastern North Pacific Ocean and identify the areas of potential concern with respect to OA in the North Pacific subpolar gyre, GoA, Bering Sea, and Amundsen Gulf. We used a combination of high-resolution models, chemical observations, particle tracking and genetic surveys of species morphometrics to determine exposure, biological responses related to sensitivity, and resilience or adaptive capacity of pteropods to OA.

Methods

Physical-Chemical Seawater Characterization

Water samples for physical and chemical parameters were collected on the GO-SHIP P16N cruise (Figure 1A), the Dyson cruise in the Bering Sea (Figure 1B) and in the Amundsen Gulf (Canadian Beaufort Sea-Marine Ecosystem Assessment, or CBS-MEA; Figure 1C). Water was collected with a CTD/rosette package equipped with sensors for temperature (°C), salinity (PSU), and oxygen (μmol kg^–1), and a rosette of Niskin bottles for collecting discrete water samples throughout the water column. Dissolved oxygen was measured using titrations. Dissolved inorganic carbon (DIC) and total alkalinity (TA) samples were collected at discrete depths throughout the water column. In the Bering Sea, TA and DIC lab analyses were conducted on the water samples where pteropods were present. DIC and TA samples were collected in 250 ml borosilicate glass bottles and preserved with mercuric chloride following Dickson et al. (2007). DIC was measured using gas extraction and colorimetric titration with photometric endpoint detection (Johnson et al., 1998). TA was determined using open-cell potentiometric titration with a full curve Gran end-point determination (Haradsson et al., 1997). Certified Reference Materials (CRMs) were analyzed with both the DIC and TA samples as an independent verification of instrument calibrations (Dickson et al., 2007; Feely et al., 2016). Precision of measurements for DIC and TA was better than 0.1 % and 0.2 %, respectively. The saturation state of seawater with respect to aragonite was calculated from the DIC and TA measurements with the CO2SYS program (Lewis and Wallace, 1998) using the dissociation constants of Lueker et al. (2000) and the boron constant of Lee et al. (2010). Based on the uncertainties in the DIC and TA measurements and the thermodynamic constants, the uncertainty in the calculated Ω_ar was approximately 0.02. The Ω_ar observed values were correlated with the pteropod shell parameters (dissolution, density, thickness).

FIGURE 1

Pteropod Sampling

The list of stations where pteropods were collected, along with the data on abundance (reported in individuals per m²), shell morphological and genetic characteristics across all the investigated stations (North Pacific subpolar gyre, GoA, Bering Sea, Amundsen Gulf) are presented in Figure 1 and Table 1. Pteropods were sampled at 12 stations ranging from 48° to 56°N in the Gulf of Alaska and North Pacific subpolar gyre on the GO-SHIP P16N cruise in May and June 2015. Pteropods were collected using 202 and 335 μm mesh Bongo nets that were obliquely trawled for 20 min across the upper 100 m of the water column at night. Pteropod samples from the Bering Sea were obtained during an oceanographic survey by NOAA Alaska Fisheries Science Center (DY17-04 and DY17-08) in spring and summer (April and August) 2017. Samples were collected using 505 μm mesh net at six different stations. All pteropods were preserved in 100% buffered ethanol. Such sampling method does not accurately represent pteropod abundances or their life stages, and thus, abundance data was not reported for the Bering Sea. During the 2018, pteropods in the Amundsen Gulf were collected using a 500 μm mesh Bongo net that was obliquely trawled for up to 15 min across five different depth strata: 0–25 m, 25–50 m, 50–100 m, 100–200 m, and > 200 m. All sampling was conducted during daylight hours and pteropods were handpicked and frozen (−80°C) for genetic and shell analyses. From the 2018 CBS-MEA sampling stations, pteropods were selected from stations at CPB, MTI and DUS transects to be analyzed for shell morphometrics, density, and thickness.

Exposure to OA in the GoA Based on the GoA-COBALT Model

To estimate seasonal, decadal, and long-term exposure to OA in the GoA, we used model output from a hindcast simulation conducted with the GoA-COBALT (Carbon, Ocean Biogeochemistry and Lower Trophic) model (Hauri et al., 2020). This model is based on the three-dimensional physical Regional Oceanic Model System (ROMS, Shchepetkin and McWilliams, 2005), the modified version of the COBALT (3PS-COBALT, Stock et al., 2014; van Oostende et al., 2018) model, and the moderately high-resolution terrestrial hydrological model by Beamer et al. (2016). Physical initial and boundary conditions were taken from the Simple Ocean Data assimilation ocean/sea ice reanalysis 3.3.1 (SODA; Carton et al., 2018). Dissolved inorganic carbon (DIC) and total alkalinity (TA) initial conditions were taken from the mapped version 2 of the Global Ocean Data Analysis Project (GLODAPv2.2016b; Lauvset et al., 2016) data set. DIC was adjusted to the corresponding hindcast year using regional anthropogenic CO₂ estimates by Carter et al. (2017). Nitrate, phosphate, oxygen, and silicate initial conditions were extracted from the World Ocean Atlas 2013 (Garcia et al., 2013). All other variables were initialized based on a climatology from a GFDL-COBALT simulation (1988– 2007), as described in Stock et al. (2014). The model was run from 1980 to 2013 and was forced with atmospheric pCO₂ observations from the Mauna Loa CO₂ record (¹last accessed 22 January 2018). Winds, surface air temperature, pressure, humidity, and radiation were forced with data from the Japanese 55-year Reanalysis (JRA55-do) 1.3 project (Tsujino et al., 2018). Modeled freshwater was brought in from numerous rivers and tidewater glaciers at a 1 km resolution and daily timestep with point-source river input via exchange of mass, momentum, and tracers through the coastal wall at all depths (Beamer et al., 2016; Danielson et al., 2020). Nutrient, DIC, and TA concentrations in the freshwater were based on available observations (Stackpoole et al., 2016, 2017). The reader is referred to Hauri et al. (2020) for more information on the model and a detailed model evaluation.

Exposure to OA in the Bering Sea

Ocean acidification (OA) exposure in the Bering Sea was estimated with model results from a 10 km horizontal resolution implementation of the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005; Danielson et al., 2011) denoted as “Bering10K”, coupled to the BEST-NPZ ecosystem model (Gibson and Spitz, 2011). Atmospheric forcing of air temperature, winds, specific humidity, rainfall, and shortwave and longwave radiation come from The Climate Forecast System Reanalysis (CFSR; Saha et al., 2010) and CVSv2 operational analysis (Saha et al., 2014). Tidal mixing and ice dynamics were also included. The BEST-NPZ ecosystem model recently underwent a substantial code revision aimed at improving model simulation of biological variables, namely primary productivity, through revisions to the ecosystem light and nutrient equations and by increasing the number of vertical layers from 10 to 30 (Kearney et al., 2020). This ecosystem model includes two phytoplankton groups (large and small), microzooplankton, large and small copepods, and euphausiids. Phytoplankton productivity consumes nutrients (nitrate, ammonium, iron) which are transferred to zooplankton biomass via grazing, or to detritus via mortality. Carbonate chemistry has also been added to model ocean acidity and elucidate underlying mechanisms of inorganic carbon cycling within different shelf regions (Pilcher et al., 2019). Boundary and initial conditions for DIC and TA were calculated from empirically derived fits with salinity using ship-based data collected during the Bering Sea Ecosystem Study (BEST) and Bering Sea Integrated Research Project (BSIERP). Freshwater runoff was modified to contain seasonally varying concentrations of DIC and TA representative of outflow from the mouth of the Yukon River. Atmospheric CO₂ varied seasonally following output from the NOAA Earth System Research Laboratory’s CarbonTracker 2016 Product (CT2016; Peters et al., 2007).

Estimating Spatial Connectivity Using Particle Tracking Model

To examine spatial connectivity in the GoA, we used hydrodynamic output from the model hindcasts described in Coyle et al. (2019). Atmospheric forcing of air temperature, winds, specific humidity, rainfall, and shortwave and longwave radiation were interpolated from the CFSR (Saha et al., 2010) and CVSv2 operational analysis (Saha et al., 2014). These hindcasts include coastal runoff (Beamer et al., 2016) as well as explicit tidal dynamics (and associated mixing). Tidally filtered weekly averaged velocities from this ROMS-based, 3-km resolution model were re-gridded to regular latitude-longitude-depth coordinates, and these values were subsequently interpolated to track particles forward in time, using a daily time step, at a constant 200 m depth.

The literature suggests that, due to prevalent currents, pteropods from the eastern North Pacific Gyre are advected (McGowan, 1963) into the GoA via northwards transport. As such, we initiated particle tracking in the offshore regions at 53°N in April for the following reasons: first, the prevalence of the northward advection that can carry pteropods from the offshore eastern Pacific subpolar gyre toward the GoA within the duration of their life cycle; second, the known offshore regional presence of pteropods and their spring spawning (e.g., Mackas and Galbraith, 2012; Doubleday and Hopcroft, 2015); third, southernmost boundary in the model domain where the initiation was feasible. The duration of particle tracking for 8 or 12 months was used to correspond to the pteropod life history, delineated by the completion of their early life stage into subadults or adults (Fabry, 1989). The particle tracking location with the shortest duration of initiation was thus used as an indicator of the early, most sensitive life stages, while the location of particles with the longest duration indicated the presence of the adults.

To examine spatial connectivity in the Bering Sea, a similar procedure was employed using output from the ROMS-based, 10-km resolution model described in Kearney et al. (2020). The physical structure and forcing of this model was identical to the version used for the Bering Sea OA exposure methods, as described above for the ROMS Bering 10K model. Regularly gridded, weekly averaged velocities from this model hindcast were interpolated to track particles backward in time using 60°N, 173°W as an ending location (on April 15, 2017) in a daily time step at a constant 50 m depth, which is an average depth distribution for L. helicina.

Morphometrics of Pteropod Morphotypes

Individuals of two different morphotypes of species Limacina helicina, recognized as L. helicina helicina and L. helicina pacifica, were selected from the samples to conduct the analyses at the morphotype level. The samples did not include morphotype characterization for individuals smaller than 0.5 mm in diameter. Individual pteropod height and width were measured with the Amscope software (Irvine, CA, United States) to subsequently determine specific morphotype. Separation of different morphotypes is possible only among the adults (McGowan, 1963; Janssen et al., 2019); thus, we only analyzed the adult life stages. Separation of the adults was first based on the presence or absence of striation (Janssen et al., 2019) and the shape of the shell (high spired shape vs flattened). The shape was reported as a height-diameter ratio (H/D) ratio (McGowan, 1963). Lower H/D ratios (flat shells) with no striation identified adult L. helicina pacifica, while higher H/D (with extended shell spire) and striation presence identified the adults of L. helicina helicina (Figure 2). When H/D ratios were plotted against the diameter of the shell, the shell proportions were fairly constant throughout the size range of the individuals collected, allowing subsequent separation of the two morphotypes (McGowan, 1963). Altogether, 98 individuals were analyzed for the morphometric analyses: 56 from the Eastern Pacific subpolar gyre, including the GoA, 17 from the Bering Sea, and 25 from the Amundsen Gulf. However, the sampling approach was not adequate to quantitatively determine spatial distribution of different morphotypes. We used both morphotypes for genetic analyses and shell dissolution analyses, while only L. helicina pacifica was used for investigating shell dissolution in combination with bulk shell density and thickness, due to insufficient sample numbers for L. helicina helicina.

FIGURE 2

Bulk Shell Density and Thickness Using Microcomputer Tomography (μXCT)

For measures of pteropod shell deposition processes, an average of 10 organisms of L. helicina pacifica per station were analyzed for bulk shell density and thickness using micro-computed tomography (μXCT; Figure 2). Only the adult animals with greater than 4 whorls were analyzed, so that shell deposition patterns could be related to a single life stage, thus excluding the early life stages that could potentially introduce life-stage bias. We differentiated two morphotypes based on their H/D ratio (McGowan, 1963; Janssen et al., 2018), as well as the absence of striation. Therefore, density and thickness measurements were morphotype-specific (only on L. helicina pacifica). A total of 83 air-dried individuals were analyzed, 41 from the Eastern subpolar gyre with the GoA, 16 from the Bering Sea and 25 from the Amundsen Gulf. After overnight drying, a μXCT ScanXmate–D160TSS105 (Comscantecno Co., Ltd., Yokohama, Japan) was used to apply X-rays to the sample stage. This system can rotate through 360 degrees, with high-resolution settings (X-ray focus spot diameter, 0.8 μm; X-ray tube voltage, 80 kV; detector array size, 1024 × 1024; 1800 projections/360°, four-times averaging). This X-ray analysis focused on the aragonite shell; therefore, the soft tissue inside the shells was not considered and was numerically removed from the 3D images. Spatial resolution of the scanned image was from 1.0 to 1.2 μm, depending on the size of the shell. Reconstruction of 3D tomography was obtained by the convolution back projection (CBP) method used the ConeCTexpress (Comscantecno Co., Ltd., Yokohama, Japan). Calculation of shell thickness and density was performed by Molcer Plus imaging software (White Rabbit Corp., Inc., Tokyo, Japan).

To evaluate the bulk density of shell, we used the CT number. The CT number is a relative number of mean grayscale value with respect to the standard material of known density. When the object is constituted by a single substance, the CT number is proportional to the bulk density. In this study, we used a limestone grain (calcite crystal, calcite standard NIST RM8544 (NBS19)) with a 300 μm diameter as the standard material, which has the same composition of the target material. This limestone particle was used on individual pteropod shells for all analysis and the calculation of the CT number. Pteropod bulk shell density was evaluated by 16-bit grayscale contrasts (65,536 gray level gradations) of a transparent image achieved by μXCT. We used the CT number of the pteropod shell which is represented by equation 1:

where μ_shell, μ_air, and μ_calciteSTD were the X-ray attenuation coefficients of the sample, calcite, and air, respectively. The bulk shell density for an entire test was evaluated with equation 2:

where n was the CT number, T_n was the total number of voxels with a specific CT number (n), and T was the total number of voxels in the whole shell. The highest density was 1000 and the lowest density was 230, being represented by the boundary between the pteropod shell and the surrounding air. The mean CT number indicated the mean density of an individual test as normalized per shell length, which we refer to as ‘normalized shell thickness’ or ‘normalized bulk shell density’ in the text below.

Bulk shell density (g per cm^–3) was calculated from the mean CT number as following:

Shell morphometrics (H/D ratio), thickness, and density data were compared with the observational data related to Ω_ar to investigate the impact of carbonate chemistry on shell-building processes over the depth of the water column inhabited by the pteropods (sensuBednaršek et al., 2014). In addition, shell density and thickness were statistically correlated to the H/D ratio across all the investigated basins using either linear or polynomial regression (Supplementary Figure 1).

Shell Dissolution

We identified the incidence of different types of shell dissolution across individuals (Table 1 and Supplementary Figure 2) collected along carbonate chemistry gradients and compared shell dissolution across different regions of the subpolar gyre, GoA and Bering Sea. While evaluating shell dissolution, it is important to have as unbiased process as possible. Thus, while we randomly selected the individuals for shell dissolution not to introduce any bias in the shell examination process, we also made sure we had a sufficient number of organisms present to determine any difference in the shell dissolution between two different morphotypes. When the number of organisms allowed for this, we aimed to equally split the analyzed subsample between two different morphotypes (in the North Pacific subpolar gyre and the GoA); however, when this was not possible due to insufficient amount of pteropod individuals (such as in the Bering Sea), then we analyzed whatever morphotypes were present. The occurrence of pteropod shell dissolution in the Amundsen Gulf, as relevant to this study, has been described by Niemi et al. (2021) and here only discussed briefly. Shell dissolution was assessed for each individual and an average value was calculated per station (Table 1 and Supplementary Figure 2). We analyzed 10-15 individuals from each of 12 stations from the subpolar gyre, including the GoA, and 6 stations from the Bering Sea.

We prepared each individual by removing the organic layer using a combination of protocols (Bednaršek et al., 2012b, 2016a). Briefly, pteropods were transferred from 100% to 70% and 50% ethanol in gradual steps, rinsed thoroughly with DI water, exposed to diluted bleach (5% sodium hypochlorite) for a treatment of approximately 30 min that was subsequently rinsed with DI water for 2-3 times to clean the shell surface, and subsequent air-dried for 12 h. Mounted samples were coated with a combination of Au-Pd for 120 s at 35A before being examined under the scanning electron microscopy (SEM).

For measures of shell dissolution severity, the presence of three categories of dissolution was assessed following Bednaršek et al. (2012b). Type I indicated a shallow, upper-prismatic level of dissolution, followed by Type II with modest dissolution protruding deeper into the cross-lamella layer, and Type III indicating the most severe dissolution protruding deeper into the crystalline layers (Supplementary Figure 2). The percentage dissolution type occurrence was assessed per each individual and correlated with the observational cruise data related to Ω_ar in the subpolar gyre, the GoA and the Bering Sea.

Genetic-Based Analyses

Sample Collection

Between five and 17 individuals per station were sampled from the Bering Sea and Gulf of Alaska cruises (Table 1). Sampling was opportunistic and depended on frequency and condition of individuals collected at each station. An additional 13 individuals were sampled from the Amundsen Gulf to act as “outlier” samples to determine species or subspecies distribution across broader oceanographic basins (Table 1). Individuals were photographed prior to DNA extraction and their morphotypes were recorded using the methods described above.

DNA Extraction, PCR, and Sequencing

DNA from each individual was extracted using the QIAGEN DNeasy tissue extraction kit (Hilden, Germany). An approximately 600 base pair region of Cytochrome c oxidase I (COI) unit was amplified using universal primers for marine invertebrates (Forward: GGTCAACAAATCATAAAGATATTGG, Reverse: TAAACTTCAGGGTGACCAAAAAATCA; Folmer et al., 1994; Geller et al., 2013; Sromek et al., 2015b). Polymerase chain reactions comprised 10–15 ng of template DNA, 4.0 mM MgCl₂ and GoTaq G2 Master Mix (Promega, United States), using an initial cycle at 98°C for 3 min, 30 cycles of 10 s at 98°C, 30 s at 55°C and 30 s at 72°C, and 5 min at 72°C. PCR products were Sanger sequenced in both directions at a commercial facility (Molecular Cloning Laboratories, San Francisco, CA, United States). Consensus haplotype sequences for each individual were created using the COI forward and reverse sequences obtained.

To determine whether there was evidence for cryptic species within the regions we sampled, we examined genetic distances between individual haplotypes and visualized the data using phylogenetic trees. The haplotype distribution across the two morphotypes was explicitly examined. To determine whether there was evidence for population structure, genetic distances were calculated within and between haplotypes at each geographic station. Finally, the sequences generated as part of the current study were compared to previously published sequences (Supplementary Tables 1, 2), to determine whether there was any evidence for species or population differentiation across the range of L. helicina.

In all analyses, sequences were aligned between individuals using the Muscle algorithm (Edgar, 2004). Pairwise distances (p-distance) between haplotypes and populations were estimated in MEGA V7.0 (Kumar et al., 2016), and standard error estimates were derived using 500 bootstrap replicates (Supplementary Tables 3, 4). Phylogenies between haplotypes were determined using a Minimum-spanning network (Bandelt et al., 1999), calculated in PopART (Leigh and Bryant, 2015).

Results

Habitats of Concern in the Subpolar Gyre and GoA as Related to OA

The P16N cruise data (from 2015) in the subpolar gyre and GoA shows that the dissolution-inducive corrosive waters (threshold for severe dissolution Ω_ar < 1.2; Bednaršek et al., 2019) are shallowest at about 53°N, 152°W (Figure 3), with the isolines of Ω_ar < 1.2 deepening both to the south in the southern subpolar region and to the north on the GoA continental shelf. The corrosive waters are characterized by “low” temperature (< 6°C), high DIC (> 2180 μmol kg−1), and Ω_ar < 1.2 (Figures 3A–F).

FIGURE 3

GoA-COBALT model output (Hauri et al., 2020) of the time series of Ω_ar at 54°N, 152°W and at 57°N, 152°W indicate that the Ω_ar shoaling into the upper 50 – 100 m of the water column occurs during winter and early spring, partially extending into summer (Figure 4). At depths greater than 100 m, most of the deeper waters across the GoA and shelf region are undersaturated with respect to aragonite (Ω_ar < 1). In the GoA model outputs from 2013, there is a west-to-east deepening of Ω_ar in the upper 100 m of the water column, with the minimum Ω_ar occurring at about 51-55°N, 152°W, in the subpolar gyre (Figures 3, 5). Accordingly, the depth of the Ω_ar saturation horizon (where Ω_ar = 1) is shallowest in the western part of the study region (< 100 m) and deepens toward the eastern side of the subpolar gyre to approximately 110-130 m (Figure 6). The highest saturation state and deepest aragonite saturation horizon occur in the summer and autumn months, and significantly lower saturation states and shallower Ω_ar horizon occur during winter and early spring months (Figure 4). Summertime subsurface Ω_ar (spatially averaged across 50 to 100 m depth, and temporally averaged across June, July, and August) is generally < 1.3 across the region, especially in the near-shore regions of the northern and eastern GoA. From Kodiak (57°N/150°W) southward, wintertime (temporally averaged across December, January, and February) subsurface Ω_ar increases to levels > 1.3 along the shelf, while off the shelf, a large area stays undersaturated with respect to Ω_ar (Figure 4A). In the central region of the subpolar gyre and the GoA, Ω_ar saturation horizon is located near 50 m throughout the year (Figure 5A), while seasonal changes in the Ω_ar saturation horizon in the eastern part of the GoA cause shorter periods of exposure to low Ω_ar conditions (Figures 5A,B).

FIGURE 4

Hauri et al. (2020).

FIGURE 5

Hauri et al. (2020).

FIGURE 6

Hauri et al. (2020).

Decadal Changes in OA Habitats in the GoA Over the 34-Year Period

In evaluating long-term change in the subpolar gyre and GoA since the 1980s, the model output suggests significant decline of Ω_ar in the upper 100 m of the water column (Figure 4). The area of undersaturated water over the 50 to 100 m depth range more than doubled between 1980 and 2013 at the location of the model output (Figure 4A). Such modelling results are consistent with large-scale field observations (Carter et al., 2017). While some subsurface waters are already undersaturated with respect to aragonite in 1980, areas of Ω_ar < 1 are now much more widespread across the GoA. Furthermore, in 1980, Ω_ar frequently reach levels < 1.2 in waters deeper than 50 m, while in more recent years, Ω_ar also frequently reaches levels < 1.2 in the upper 50 m of the water column in both on- and off-shelf areas on a seasonal basis (Figure 6). In particular, the subpolar gyre is significantly lower in Ω_ar in the winter months. While Ω_ar conditions in the subpolar gyre in the 1980s ranged from 0.8 < Ω_ar < 1.2, current ranges are from 0.5 < Ω_ar < 0.8. Overall, the center of the subpolar gyre is currently the most vulnerable to the effects of increasing OA (Hauri et al., 2021). Although the eastern part of the GoA has seen a similar magnitude of change, the current OA conditions are less severe.

Current Habitats of Concern With Respect to OA in the Bering Sea

The Bering Sea 10K model output suggests that in the subsurface water between 50 and 100 m Ω_ar is < 1.3 across nearly all of the Bering Sea shelf, and Ω_ar < 1 for large portions of the northwestern shelf for both winter (December-February) and summer (June-August) months (Figure 7). In Figure 7, regions with a depth < 50 m, which are representative of the inner-shelf domain, are excluded. Subsurface Ω_ar values are generally lower in summer months compared to winter months, except for a relatively small region of buffered waters within Bristol Bay. The most corrosive waters are located in the northwestern Bering Sea. Panels A and B in Figure 8 illustrate a depth profile of model output for Ω_ar along the midnorth (MN) shelf transect (denoted by the magenta line in Figure 7). Phytoplankton productivity in surface waters and respiration at depth generate a very strong vertical gradient in Ω_ar during summer, compared to a relatively weak vertical gradient in winter. Surface Ω_ar values are greater than 2 during the summer, but less than 1.3 during the winter. Panel C in Figure 8 shows the modeled vertical profile of Ω_ar at an outer shelf location over the entire 2003–2012 period. This figure illustrates the seasonal cycle of low winter surface values, followed by higher surface Ω_ar values in spring-summer, and a subsequent subsurface drawdown in autumn due to respiration. Such intense seasonal changes in the spatial and vertical Ω_ar conditions define the exposure time of less favorable conditions (Ω_ar < 1.3) on a seasonal basis. The subsurface values of Ω_ar tend to decrease by ∼0.1 over the 2003-2012 timeframe, coinciding with a shift in physical variability from a warm temperature period to a cold temperature period. This temperature shift in the colder years due to increased vertical mixing and nutrient resupply in the euphotic layer (Pilcher et al., 2019) generates an increase in primary productivity and anthropogenic carbon uptake, with subsequent increase in respiration that exacerbates the expected OA trend over this timeframe. The surface trends of 0.025–0.04 units/year are significant for the inner and middle shelf domains (Pilcher et al., 2019).

FIGURE 7

FIGURE 8

Current Conditions in the Amundsen Gulf

Generally, in the Canadian sector of the Amundsen Sea, the supersaturated conditions occurred in surface waters, while undersaturated conditions (Ω_ar < 1) occurred across the Gulf in the Pacific halocline layer (50–200 m). Some smaller-scale variation in the Ω_ar is possible, depending on the transects sampled (Figures 9A–C), with such conditions described in more detail by Niemi et al. (2021). In the areas where pteropods were collected, Ω_ar in the Pacific layer ranged from 0.76 to 1.3, and averaged around 0.88, indicating that the pteropods’ main habitat is largely corrosive in summer (Figure 9).

FIGURE 9

Sensitivity-Related Biological Parameters

Morphotype Distribution

Pteropods were present across all of our study regions in high abundances, especially from 55 °N (Table 1). Based on the presence of striations, as well as H/D ratios, two adult morphotypes were detected: L. helicina helicina with an inflated spire and striated whorls (Figure 2D), and L. helicina pacifica with a depressed flat spire and no striation (Figure 2F). The H/D ratios plotted against the diameter of the shell showed a relatively constant ratio, regardless of the size range of the adults (Supplementary Figure 1A). The ranges of H/D values for the two morphotypes were similar to those previously reported by McGowan (1963). The H/D range for L. h. pacifica was between 0.8 to 1, while H/D range for L. h. helicina was higher than 1. All morphotypes were found in all three basins (Supplementary Figure 1B).

Shell Thickness and Density Related to the Shell Morphometrics

The distribution of H/D ratios in two morphotypes was also related to a specific range of normalized shell thickness and density (Supplementary Figures 1C,D). Normalized bulk shell thickness showed an inverted bell curve shape relationship with the H/D ratio; its decrease was correlated with declining H/D ratios to a value of approximately 0.9 to 1. Below this ratio, normalized shell thickness started to increase (polynomial fit with R² = 0.17, p-value = 8.622E-05). These results indicate that the organisms with lower spire (lower H/D range) built thicker shells compared to those on the upper range that built thinner shells (Supplementary Figure 1C). The density was not dependent on the H/D ratio, with similar densities across the range (R² = 0.052; p = 0.0155; Supplementary Figure 1D).

Environmental Control on Shell Building Processes

The impact of Ω_ar on the processes of pteropod morphometrics and shell characteristics (H/D ratio, normalized bulk shell density and thickness) across all investigated ocean regions (GoA, Bering, Amundsen) was examined in one morphotype, L. h. pacifica only. The H/D ratio of adult individuals was significantly correlated with Ω_ar. Linear regressions showed the best fit of H/D ratio to be correlated with the Ω_ar, at 50 m (Ω_ar,50) representing an average depth of daily pteropod vertical migration (Figure 10A; R² = 0.6966, p = 8.074E-16). Linking shell morphometrics with the Ω_ar, adult individuals with higher H/D ratios (extended shell with higher spire) were observed at higher Ω_ar, while the organisms with lower H/D (flatter shell with lower spire) were present at lower Ω_ar. We fitted normalized bulk shell thickness against Ω_ar averaged over 100m depth (Ω_ar,avg100). The best linear fit shows the shell thickness to significantly increase with a declining Ω_ar,avg100 (Figure 10B; R² = 0.31, p-value = 0.028). At higher values of Ω_ar, adult pteropods built shells with higher H/D ratios and lower shell thicknesses (Figures 10A,B; Supplementary Figure 1). We also fitted shell thickness against the Ω_ar at 100 m; the best fit found was polynomial one, where shell thickness increased up to a Ω_ar of 1.1, upon which it started to decrease (R² = 0.623; p-value = 2.581E-07; Figure 10C). This value could be considered an inflection point (threshold), below which the shell building processes start to decline. Such patterns were observable when combining the data of L. helicina pacifica across investigated ocean regions. However, on the regional basis, the lack of sufficient data or limited Ω_ar,avg100 condition gradients prevented us from making conclusions for each region separately, especially for the Bering Sea region with its narrow range of Ω_ar,avg100 conditions. Normalized bulk shell density was not significantly correlated with Ω_ar,avg100 (R² = 0.027, p-value = 0.273; Figure 10D). These results imply that only the processes influencing shell thickness, and not density, were impacted by the carbonate chemistry conditions.

FIGURE 10

Pteropod Shell Dissolution Across the Regions

High resolution observations of pteropod shell surface under SEM demonstrated a range of dissolution values, depending on observed Ω_ar conditions at the specific sample locations (Figure 2). The samples collected during spring-summer in the subpolar gyre and the GoA were characterized by a notable difference in the shell dissolution, with the stations across the central and western side of the study region (stations 178-191, Figure 11) showing the greatest severity of dissolution (Table 1). The areas of dissolution spatially coincide with the regions characterized by a shallower Ω_ar depth horizon, depicted in both the observational data and model outputs (Figure 5). In comparison, shell dissolution declined toward the east, with complete disappearance of any dissolution patterns east of 139°W. In addition, the more inshore locations (Station 187) had less severe dissolution compared to the offshore samples. This result again coincides with the GoA-COBALT model output of higher Ω_ar onshore and lower Ω_ar and shallower Ω_ar saturation horizon offshore (Figures 5, 11).

FIGURE 11

Out of the six collected samples in the Bering Sea during spring-summer, various dissolution ranges were observable (Table 1), mostly characterized by Type I and II dissolution, with only station 24 having more severe (Type III) dissolution present. Because there was insufficient data from the cruise to interpolate across the stations, we could not correlate shell dissolution with the oceanographic gradients as in the GoA, so we used the model output instead (Figures 7, 8). The stations with the most severe dissolution were distributed over the central part of the eastern Bering Sea shelf. This coincides with the model outputs of the observed shallow Ω_ar horizon and Ω_ar < 1 around 50 m depth during this period, which is also in the Bering 10k model output. Intact organisms were present in the southern part of the Bering Sea and around the Aleutian Islands, with the Ω_ar horizon positioned at depth ranges below 100 m.

Shell dissolution in the Amundsen Gulf was described in detail by Niemi et al. (2021). In short, shell dissolution was widespread, occurring at coastal and offshore locations as well as in small embayments. In an initial assessment of 134 samples of L. helicina, greater than 85% of individuals exhibited severe dissolution. The majority of shell damage occurred on the first whorl, indicating that damage likely occurred in early life stages during the spring period (May/June) when sea ice may still be present in the Amundsen Gulf. Overall, we found no observable difference in shell dissolution between two different morphotypes across any of the high latitude regions, indicating that shell dissolution is not a trait that would be morphotype-specific.

Biological Parameters Related to Resilience and Adaptive Capacity

Genetic Analyses

A total of 176 individuals were sequenced at cytochrome c oxidase I (COI) across 21 stations (Supplementary Tables 1, 2). Of these, 130 represented the L. helicina helicina phenotype and 38 represented the L. helicina pacifica phenotype. The phenotype of eight sequences were not determined due to sample deterioration. After alignment, 544bp were shared across all individuals, representing 20 unique COI haplotypes across the entire dataset.

Haplotype diversity was low (Supplementary Table 3): overall nucleotide diversity was 0.005 (SE 0.001). Pairwise distances between individual haplotypes varied between 0.002 and 0.009. These values were lower than accepted for species definition at the CO1 locus (Hebert et al., 2003). The samples were dominated by three haplotypes (Figure 12A): Haplotype 1 was found in 118 individuals, Haplotype 2 in 34 individuals, and Haplotype 3 in 5 individuals. The remaining haplotypes were found in two individuals (Haplotype 4 and 5) or in one individual only (Haplotypes 6-20).

FIGURE 12

The minimum spanning network described a star phylogeny for CO1 haplotypes (Figure 12A); that is, several short branches connected by a common haplotype. Such phylogenies are often interpreted to represent an evolutionarily recent founder event, followed by population expansion. There was no evidence for genetic differentiation between L. helicina morphotypes (Figure 12A). Both morphotypes were represented within three most frequent haplotypes.

Visualization of the frequency distributions of the three most common haplotypes (1-3) suggests that sequence diversity decreases in the northernmost populations, particularly in the Bering Sea and Amundsen Gulf (Figure 12B). However, formal analysis of haplotype diversity between stations sampled revealed no significant differences between populations (Supplementary Table 4). Divergence of haplotype frequencies between populations varied between 0.001 and 0.002, with standard errors falling in the same order of magnitude. The power of this analysis may be affected by sample sizes, which varied between 5 and 17 individuals per station. However, given the high frequency of the most common haplotype, as well as the close relationship between these haplotypes, large sample sizes per station would be required to detect putative differences.

A total of 259 published COI sequences (Supplementary Table 2) were combined with the 179 sequences derived from this study and were grouped into nine broad geographic regions: Amundsen Gulf, Bering Sea, subpolar gyre, Greenland, Kara Sea, Hudson Bay (Northern Canada), Western North Pacific, Svalbard and White Sea. Global haplotype divergence was small, with most haplotypes differing from each other by a few mutations (Supplementary Figure 3). However, two previously described haplogroups (Abyzova et al., 2018; Shimizu et al., 2018) described the data, with evidence of a unique haplogroup limited to Svalbard and the Kara Sea, and dominance of the second haplogroup in the Bering Sea, Amundsen Gulf and North Pacific Seas. While the haplogroups differ by only a few mutations, these data suggests separate colonization and divergence events and broad-scale population structure. However, pairwise tests for divergence (Supplementary Table 5) revealed small pairwise distances between geographic regions.

Spatial Connectivity Based on Particle Tracking Model Outputs

Particle tracking models were used to investigate pteropod spatial connectivity over the subpolar gyre in the GoA. These were based on the assumption of a northward juvenile pteropod advection from the eastern North Pacific and life duration of approximately eight months to a year, during which pteropods would reach adulthood. The results indicate that pteropods were dispersed throughout the entire area during this period, irrespective of location origin (Figure 13), with high abundances present (Table 1). On the other hand, shorter durations of particle tracking indicated the locations of the early life stages, which are dependent on the location of their spawning.

FIGURE 13

Coyle et al. (2019) and Kearney et al. (2020).

Northward advection dominates in the first three months after the typical spawning period in late spring. In the following 4–8 months, flow of the GoA waters occurs westward along the Aleutian Islands (Figure 13A). If initiation is assumed to occur in the offshore subpolar gyre, the early life stages can thus reach the GoA in the summer, while the subadults and adults would be distributed in the GoA’s coastal and shelf region and advected westward in the winter. The dispersal of the particles over long distances demonstrates large-scale biogeographic distribution with no apparent species boundaries across the entire area.

Extensive spatial connectivity was also observed in the Bering Sea (Figure 13B), where we tracked particles backward in time for a full year from a mid-shelf location on the northwestern shelf. This ending location was chosen based on observational results. Particles reaching this destination all originated from the south, primarily from the Gulf of Alaska (passing through Unimak Pass), but also from the deep basin to the north of the Western Aleutian Islands. From either starting location the particles travel along the shelf break to the north of the Aleutian Islands during an approximately 8-month period, subsequent to their initial release in April. In general a smaller spread of particles was observed for the Bering Sea than was the case in the GOA. Note these results are not strictly comparable, as we tracked particles forward in time in the GOA case, and backwards in the Bering Sea case. It is further possible that more dispersion would have been obtained with a different ending location in the Bering Sea. Fundamentally, however, the smaller spread of the Bering Sea tracks reflects the presence of vigorous 200-km scale eddies in the deep GOA, as compared to the relatively quiescent conditions on the shallow Bering Sea shelf.

Discussion

This interdisciplinary study provides an integrated OA risk assessment for high-latitudinal pteropod populations of the subpolar gyre, GoA, Bering Sea and Amundsen Gulf. Here, we synthetize multi-faceted components of the risk assessment, as defined by OA exposure, species and morphotype sensitivity, and resilience through genetic structure and spatial connectivity.

Model output comparisons between the two regions (the GoA and the Bering Sea) revealed the differences in the magnitude of Ω_ar conditions, both in the context of the vertical gradients and seasonal cycles related to Ω_ar. Overall, the model outputs as well as data based on P16N cruise observations (Figure 1), suggest that Ω_ar is lower in the North Pacific subpolar gyre, while the habitats in the eastern part of the GoA seem to be less OA-compressed, providing more favorable conditions for pteropod populations. The Bering Sea, in comparison, has somewhat lower values of Ω_ar and a shallower Ω_ar horizon than the subpolar gyre/GoA. This observation is consistent with decreasing Ω_ar towards the poles in North-Pacific-Arctic waters (Mathis et al., 2015), indicating that this region may have experienced more severe biological impacts corresponding with the changes, especially in the last two to three decades (Carter et al., 2017). Seasonality is an important parameter to consider, because it defines the duration of exposure to less favorable conditions. In the GoA, the subsurface conditions are less favorable in the summer, while the subsurface conditions in the Bering Sea are less favorable in the winter (Figures 5–8). As illustrated by the Bering Sea model, vertical habitat suitability is reduced in the summertime, showing persistent subsurface Ω_ar < 1.2 conditions despite high surface Ω_ar conditions that are buffered by substantial primary productivity. There is a clear exception in the subpolar gyre, where unfavorable conditions appear to persist throughout the year. Still, while surface waters of the subarctic Pacific are projected to become Ω_ar < 1 (Bindoff et al., 2019), global model simulations suggest that surface waters of the Bering Sea will remain supersaturated even under the RCP 8.5 scenario. Hence, the upper surface waters represent suitable habitat for sensitive pelagic calcifiers. However, it remains uncertain if the vertical migrators, such as pteropods, could only occupy surface waters given their ecological dependence on deeper depths (Lalli and Gilmer, 1989).

With respect to the trend observed over the last 34 years, the GoA-COBALT model output of the subpolar gyre and the GoA suggests a declining trend of Ω_ar and shoaling of the Ω_ar saturation horizon. This indicates that pteropod populations have been experiencing an increase in exposure to lower Ω_ar in their habitats over decadal time scales. While we do not have biological data over the same time period as the model output, we assume that it is very likely that the severity of pteropod shell dissolution described here is a consequence of rapid, multiple human-related impacts in the region, and not a natural change in the region (Mathis et al., 2015; Carter et al., 2017; Hauri et al., 2021). The lack of a hindcast simulation in the Bering Sea does not allow such biological risk assessment over a longer time period.

Insights into the seasonal carbonate chemistry variability is important to infer the impacts of exposure within the L. helicina life history. The life cycle of northern high-latitude pteropods is regionally specific, lasting between 1 and 2 years (Kobayashi, 1974; Fabry, 1989; Gannefors et al., 2005), and is characterized by late spring-summer spawning, with subsequent presence of subadults and adults during the winter periods of the corresponding years. Successful recruitment depends mostly on survivorship during two key periods: the early life stages that have higher sensitivity to OA, and the adult stages that are responsible for successful spawning and recruitment (Bednaršek et al., 2016b). Model outputs during the spring-summer period can be used to characterize exposure of the early life stages (larval and juveniles) to OA, while the model outputs integrated over the upper 50–100 m of the water column for the winter months can serve as an approximation of the adult exposure (Bednaršek et al., 2019). The summer conditions in the GoA are more favorable, thus early life stages will be less impacted relative to the adults that experience lower Ω_ar wintertime exposure. In contrast, the Bering Sea has opposite seasonality, characterized by the summer reduction of the suitable Ω_ar conditions, with adults being less impacted compared to the early life stages. It is also important to note that regionally specific conditions, like the year-round unfavorable conditions in the subpolar gyre, might impact both early and adult life stages, delineating this region of the highest OA concern.

Observations of adult pteropod shell dissolution during spring-summer surveys in both basins spatially correlate with the spring spatial carbonate chemistry patterns, as delineated by the P16N cruise observations and biogeochemical models in the GoA and subpolar gyre (Figures 1, 11). In the GoA, there was more severe dissolution in the western part of the study area, coinciding with the highest pteropod abundances, and less severe dissolution in the onshore waters. In the Bering Sea, severe dissolution occurred in the north-central part, where the model outputs projected lower Ω_ar compared to the higher Ω_ar in the south. The GoA is more likely to experience severe spring and summertime dissolution of pteropod shells compared to the same period in the Bering Sea. Our data show that thinner shells, along with dissolution uniformly distributed around the shell surface, are evidence of an impact related to lower Ω_ar cumulative exposure (Bednaršek et al., 2016a). This is opposite to the conclusion from Peck et al. (2016); namely that dissolution only occurs if the periostracum is damaged or if imposed through predation pressure. More severe dissolution in the Amundsen Gulf likely occurs because of the dominance of undersaturated Pacific waters within the pteropod habitat, where prolonged low Ω_ar exposure could represent higher risks, characteristic of the polar habitats of the Northern Hemisphere. Given that our genetic analyses show that the Amundsen Gulf pteropods share the same species distribution as the GoA and the Bering Sea individuals, there is no reason to believe that another more resilient species or subspecies will replace the current populations (Niemi et al., 2021).

There are additional parameters that can impact shell processes (Bednaršek et al., 2018; Oakes and Sessa, 2020), including regional parameters that can reverse the shell patterns observed in this study (Mekkes et al., 2021), making it important to compare the differences always among the same morphotypes. Seasonal patterns of high primary production in the Bering Sea could potentially offset some of the negative effects (sensuBednaršek et al., 2018), with temperature and oxygen also playing a role in the shell building processes (Mekkes et al., 2021). Both the GoA and the Bering Sea pteropods have substantially less severe dissolution compared to that reported in the Arctic regions, such as the Amundsen Gulf. Limited extent of shell repair in the form of crystalline regrowth was reported in the presence of abundant food availability only in a few adult organisms (Niemi et al., 2021). Against a general pattern of thinner shells with only a few individuals displaying limited shell repair, this data shows high risk for pteropod population. Our findings therefore do not support the observation of pteropod resilience to OA by Peck et al. (2018).

Sensitivity and adaptation capacity were evaluated by examining the shell morphometric characteristics. Flatter and thicker (lower H/D) individuals were found in habitats associated with more severe Ω_ar, while high-spire, less thick (high H/D) individuals were present at higher Ω_ar. This distribution likely indicates that pteropods can build thicker shells and display some plastic responses, thus indicating a phenotypic response in the less favorable conditions that protects them against intense shell dissolution. A divergent range of plastic responses in the two distinct phenotypic morphotypes could be differential acclimation to Ω_ar exposure. The extended height over the same diameter results in higher surface shell area of L. helicina helicina, pointing towards the fact that a larger extent of shell surface can be exposed to the surrounding waters compared to L. helicina pacifica with lower surface area. As such, we postulate that different morphotypes might develop in the native habitats characterized by the intensity of Ω_ar exposure. By inhabiting more intense OA conditions, L. helicina pacifica may form flatter and thicker shells to protect the shell against more severe exposure to lower Ω_ar. Flatter shells with lower surface exposure is intended to reduce the exposure to low Ω_ar, while thicker shells can potentially tolerate more extensive dissolution better at lower Ω_ar. On the other hand, L. helicina helicina inhabits the habitat with less severe, low Ω_ar exposure and does not need such protecting mechanisms; greater shell surface and thinner shells can be sustained at less intense OA condition.

Pteropod morphotype distribution data over much larger spatial scales in the North Pacific clearly demonstrates the eastward occurrence of L. helicina pacifica, while L. helicina helicina is dominant westwards (McGowan, 1963). In combination with the carbonate chemistry data over the examined area that shows a progressing shoaling of carbonate chemistry eastwards (Carter et al., 2017; Cai et al., 2020), there is a spatial pattern of L. helicina helicina inhabiting the westward region with higher Ω_ar and L. helicina pacifica present in the lower Ω_ar eastward region, to which they seem to be more adapted to via their morphometric characteristics. However, given the presences of both morphotypes across the stations of the investigated area of the subpolar gyre and the GoA, our study is likely not spatially extensive enough to be able to demonstrate differential morphotype occurrence over Ω_ar scales.

In terms of shell building processes, pteropods continue building their shells even under lower Ω_ar conditions (sensuComeau et al., 2010) but at reduced growth rates, presumably due to the trade-offs in the compensatory processes. However, it seems that the conditions of around 1.1 < Ω_ar > 1.3 (Figures 10B,C) could be delineated as a threshold below which the shell building process becomes energetically more expensive, and thus, shell processes start to decline. A similar Ω_ar biomineralization threshold was previously reported by Bednaršek et al. (2019). Interestingly, shell density seems to be an independent process not driven by Ω_ar, likely pointing towards the fact that building a shell with consistent density is an important process for the organisms and is thus more tightly regulated than processes governing shell thickness. While the differences across the eastern-western area in aragonite saturation state were substantially different, we note minimal differences in temperature across the same area, with only 1-2°C difference (Figures 3A,B), suggesting a relatively homogenous habitat temperature wise (sensuMcGowan, 1963). As Ω_ar and temperature are highly colinear in both basins, it makes it impossible to differentiate which parameter is a major driver behind the biological responses. In addition, the lack of data on other potentially important parameters, such as food availability (sensuBednaršek et al., 2018) is currently not allowing for a more accurate interpretation and prediction of the drivers and mechanisms behind the biomineralization processes and population-level outcomes under future OA intensification.

Beside a direct impact of OA on reproduction (sensuManno et al., 2016), there could be indirect implications of energetic expenditure related to reproductive outputs. The assumption is that higher Ω_ar conditions would allow individuals to potentially invest more energy into growth, resulting in larger shells with greater tissue mass. Building shells under low Ω_ar conditions (below 1.1) would require greater energetic investment, which may be indirectly reflected in smaller-sized organisms. Since the reproductive efforts of the individuals are inherently linked to tissue weight (Lalli and Gilmer, 1989), smaller-sized organisms with lower gamete production might have a lower reproductive effort under future OA intensification.

Examining the parameters related to the adaptation capacity (genetic structure, cryptic species, and spatial connectivity) is also important to potentially determine the morphotype-level risks on pteropod population level. For example, OA-induced exposure would not necessarily result in local vulnerability if individuals were replaceable or if spatial connectivity may bring new flux of organism from less impacted areas, while the presence of genetically differentiated subspecies with differential response to OA might buffer changes in distribution and abundance. Despite morphological differences, genetic analyses reveal that there was no basis for differentiation between these two forms at the species level, because both share mitochondrial DNA haplotypes. In addition, particle tracking modeling indicates that the organisms are well mixed over the investigated area, supporting genetic uniformity of the investigated pteropods. Therefore, we suggest that the differences in morphotypes might reflect phenotypic plasticity related to the differential exposure to prevailing OA-related environmental conditions, but likely not limited to Ω_ar only. Alternatively, such differences may arise through fine- scale population structure and local adaptation; this explanation cannot be ascertained with our available data. We conclude that there is no evidence for cryptic species or subspecies across the geographic ranges sampled, based on tests for divergence at the mitochondrial DNA cytochrome oxidase I locus. It is unlikely that hidden species diversity may account for differential morphotype responses to low Ω_ar. It is possible that there is within-species diversity that should be examined using nuclear DNA markers. Our results reveal a possible reduction in haplotype diversity in northern populations, which may point toward population differentiation. However, this result was not significant. On a species-wide range, the previously reported differentiation between unique haplogroups in the Arctic Ocean north of mainland Europe (Kara Sea and Svalbard), the Western North Pacific (Shimizu et al., 2018), and the Okhotsk Sea (Shimizu et al., 2021) was supported here by the addition of further samples from the subpolar gyre, the Bering Sea, and Amundsen Gulf. While the differences between the two haplogroups were small, the results suggest population structure at the species-wide scale, explained by separate colonization events with subsequent divergence over broad geographic ranges. These results are intriguing, given that the single molecular marker used has relatively little power to differentiate between populations.

Investigations based on particle analyses and prevailing advective patterns suggest extensive spatial connectivity. These analyses further support the result of a single cosmopolitan species of L. helicina, revealed by the mitochondrial DNA analyses. Given the time trajectories of the particles modeled here, the results suggest a continuous process of repopulation and intermixing of pteropods in the GoA region in the spatial context of the North Pacific subpolar gyre. High connectivity might alleviate localized negative impacts, especially if different morphotypes are constantly repopulating this region, therefore providing resilience. However, poor survival at different life history stages may affect recruitment throughout the region. For example, the early life stages encounter exposure to low Ω_ar conditions between 53 and 56°N through northward advection. It is also unknown to what degree adaptive capacity at the population level may buffer recruitment variation. Overall, high risk for pteropod population persistence with intensification of OA has to be assumed because of their known vulnerability, along with low evidence for species diversity and the progressive severity of OA exposure, despite their wide connectivity across the considered habitats and our insufficient knowledge on their adaptive capacity in the high latitudes.

Corresponding with pteropod high risk related to OA in the high latitudinal habitats are also high-risk repercussion for the commercially important fisheries and high trophic levels marine vertebrates. Considering the magnitude of OA change in this region coupled with a high pteropod biomass (Bednaršek et al., 2012a; Doubleday and Hopcroft, 2015), negative pteropod population implications could ultimately have an impact for the fisheries dependent on pteropods as a food source, most explicitly for the various salmon species, cod and char, during the fall time in the absence of the other food sources essential for successful fish development. Moreover, given the changes in pteropod shell building processes, i.e. reduced calcification and shell thickness and lower H/D ratio, along with increased dissolution, these high latitudinal regions’ with one of the highest pteropod abundances on a global scale (Table 1; Bednaršek et al., 2012a) might also experience potentially significant biogeochemical effects related to the biological pump, resulting in the lower carbonate export fluxes over time (Carter et al., 2021; Lee and Feely, 2021; Sulpis et al., 2021). The results of this study suggest that an important component of the alkalinity change in the North Pacific may be the result of changes in pteropod production and dissolution. Moving forward into the future, the ecological and biogeochemical implications of these long-term changes will depend on multiple interacting processes, including the morphotype distribution, net shell calcification and dissolution, changes in the morphometrics, their abundance and phenology.

Given the dominance of pteropods in the subpolar gyre and the GoA, long-term observations should be strategically positioned along the W-E gradients where the differences in OA are most pronounced. Understanding of biological vulnerability on the regional scales provides an opportunity to prioritize monitoring to manage, and thus potentially reduce, the realized impact on the fisheries and ecosystems in the high-latitudinal ecosystems under predicted climate change scenarios.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  AbyzovaG. A.NikitinM. A.PopovaO. V.PasternakA. F. (2018). Genetic population structure of the pelagic mollusk Limacina helicina in the Kara Sea.PeerJ6:e5709.

  • Google Scholar
  • View reference in article

• 2

  ArmstrongJ. L.BoldtJ. L.CrossA. D.MossJ. H.DavisN. D.MyersK. W.et al (2005). Distribution, size, and interannual, seasonal and diel food habits of northern Gulf of Alaska juvenile pink salmon, Oncorhynchus gorbuscha.Deep Sea Res. Part II52247–265. 10.1016/j.dsr2.2004.09.019

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AuburnM. E.IgnellS. E. (2000). Food habits of juvenile salmon in the Gulf of Alaska July–August 1996.N. Pac. Anadr. Fish. Comm. Bull.289–97.

  • Google Scholar
  • View reference in article

• 4

  BandeltH. J.ForsterP.RöhlA. (1999). Median-joining networks for inferring intraspecific phylogenies.Mol. Biol. Evol.1637–48.

  • Google Scholar
  • View reference in article

• 5

  BeamerJ. P.HillD. F.ArendtA.ListonG. E. (2016). High-resolution modeling of coastal freshwater discharge and glacier mass balance in the Gulf of Alaska watershed.Water Resour. Res.523888–3909. 10.1002/2015wr018457

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BeauchampD. A.CrossA. D.ArmstrongJ. L.MyersK. W.MossJ. H.BoldtJ. L.et al (2007). Bioenergetic responses by Pacific salmon to climate and ecosystem variation.N. Pac. Anadr. Fish. Comm. Bull.4257–269.

  • Google Scholar

• 7

  BednaršekN.MožinaJ.VogtM.O’BrienC.TarlingG. A. (2012a). The global distribution of pteropods and their contribution to carbonate and carbon biomass in the modern ocean.Earth Syst. Sci. Data4167–186. 10.5194/essd-4-167-2012

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BednaršekN.TarlingG. A.BakkerD. C.FieldingS.CohenA.KuzirianA.et al (2012b). Description and quantification of pteropod shell dissolution: a sensitive bioindicator of ocean acidification.Glob. Chang. Biol.182378–2388. 10.1111/j.1365-2486.2012.02668.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BednaršekN.JohnsonJ.FeelyR. A. (2016a). Comment on Peck et al: vulnerability of pteropod (Limacina helicina) to ocean acidification: shell dissolution occurs despite an intact organic layer.DSRII12753–56. 10.1016/j.dsr2.2016.03.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BednaršekN.FeelyR. A.ReumJ. C. P.PetersonB.MenkelJ.AlinS. R.et al (2014). Limacina helicina shell dissolution as an indicator of declining habitat suitability owing to ocean acidification in the California Current Ecosystem.Proc. R. Soc. B Biol. Sci.281:20140123. 10.1098/rspb.2014.0123

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BednaršekN.HarveyC. J.KaplanI. C.FeelyR. A.MožinaJ. (2016b). Pteropods on the edge: cumulative effects of ocean acidification, warming, and deoxygenation.Prog. Oceanogr.1451–24. 10.1016/j.pocean.2016.04.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  BednaršekN.FeelyR. A.TolimieriN.HermannA. J.SiedleckiS. A.WaldbusserG. G.et al (2017). Exposure history determines pteropod vulnerability to ocean acidification along the US West Coast.Sci. Rep.7:4526. 10.1038/s41598-017-03934-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  BednaršekN.FeelyR. A.BeckM. W.GlippaO.KanervaM.Engström-ÖstJ. (2018). El Niño-related thermal stress coupled with upwelling-related ocean acidification negatively impacts cellular to population-level responses in pteropods along the California current system with implications for increased bioenergetic costs.Front. Mar. Sci.5:486. 10.3389/fmars.2018.00486

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  BednaršekN.FeelyR. A.HowesE. L.HuntB. P. V.KessouriF.LeónP.et al (2019). systematic review and meta-analysis toward synthesis of thresholds of ocean acidification impacts on calcifying pteropods and interactions with warming.Front. Mar. Sci.6:227. 10.3389/fmars.2019.00227

  • CrossRef
  • Google Scholar

• 15

  BindoffN. L.CheungW. W. L.KairoJ. G.ArísteguiJ.GuinderV. A.HallbergR.et al (2019). “Changing ocean, marine ecosystems, and dependent communities,” in IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edsPörtnerH.-O.RobertsD. C.Masson-DelmotteV.ZhaiP.TignorM.PoloczanskaE. (Geneva: IPCC).

  • Google Scholar
  • View reference in article

• 16

  BrownZ. W.van DijkenG. L.ArrigoK. R. (2011). A reassessment of primary production and environmental change in the Bering Sea.J. Geophys. Res. Oceans11610.1029/2010JC006766

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  BurridgeA. K.GoetzeE.RaesN.HuismanJ.PeijnenburgK. T. C. A. (2015). Global biogeography and evolution of Cuvierina pteropods.BMC Evol. Biol.15:39. 10.1186/s12862-015-0310-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  CaiW.-J.XuY.-Y.FeelyR. A.WanninkhofR.JönssonB.AlinS. R.et al (2020). Controls on surface water carbonate chemistry along North American ocean margins.Nat. Commun.11:2691. 10.1038/s41467-020-16530-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  CarterB. R.FeelyR. A.WanninkhofR.KouketsuS.SonnerupR. E.PardoP. C.et al (2019). Pacific anthropogenic carbon between 1991 and 2017.Global Biogeochem. Cycles33597–617. 10.1029/2018gb006154

  • CrossRef
  • Google Scholar

• 20

  CarterB. R.FeelyR. A.LauvsetS. K.OlsenA.DeVriesT.SonnerupR. (2021). Preformed properties for marine organic matter and carbonate mineral cycling quantification.Global Biogeochem. Cycles35:e2020GB006623. 10.1029/2020GB006623

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  CarterB. R.FeelyR. A.MeckingS.CrossJ. N.MacdonaldA. M.SiedleckiS. A.et al (2017). Two decades of Pacific anthropogenic carbon storage and ocean acidification along Global Ocean Ship-based Hydrographic Investigations Program sections P16 and P02.Global Biogeochem. Cycles31306–327. 10.1002/2016GB005485

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  CartonJ. A.ChepurinG. A.ChenL. (2018). SODA3: a new ocean climate reanalysis.J. Clim.316967–6983. 10.1175/JCLI-D-18-0149.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  ComeauS.GorskyG.JeffreeR.TeyssiéJ.-L.GattusoJ.-P. (2009). Impact of ocean acidification on a key Arctic pelagic mollusc (Limacina helicina).Biogeosciences61877–1882. 10.5194/bg-6-1877-2009

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  ComeauS.GorskyG.AlliouaneS.GattusoJ. P. (2010). Larvae of the pteropod Cavolinia inflexa exposed to aragonite undersaturation are viable but shell-less.Mar. Biol.1572341–2345. 10.1007/s00227-010-1493-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  CoyleK. O.HermannA. J.HopcroftR. R. (2019). Modeled spatial-temporal distribution of productivity, chlorophyll, iron and nitrate on the northern Gulf of Alaska shelf relative to field observations.Deep Sea Res. Part II165163–191. 10.1016/j.dsr2.2019.05.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  CrossJ. N.MathisJ. T.BatesN. R.ByrneR. H. (2013). Conservative and non-conservative variations of total alkalinity on the southeastern Bering Sea shelf.Mar. Chem.154100–112. 10.1016/j.marchem.2013.05.012

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  DalyE. A.MossJ. H.FergussonE.DebenhamC. (2019). Feeding ecology of salmon in eastern and central Gulf of Alaska.Deep Sea Res. Part II165329–339. 10.1016/j.dsr2.2019.06.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  DanielsonS.CurchitserE.HedstromK.WeingartnerT.StabenoP. (2011). On ocean and sea ice modes of variability in the Bering Sea.J. Geophys. Res. Oceans116:C12034. 10.1029/2011JC007389

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  DanielsonS. L.HillD. F.HedstromK. S.BeamerJ.CurchitserE. (2020). Demonstrating a high-resolution Gulf of Alaska ocean circulation model forced across the coastal interface by high-resolution terrestrial hydrological models.J. Geophys. Res. Oceans125:e2019JC015724. 10.1029/2019JC015724

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  DarnisG.BarberD. G.FortierL. (2008). Sea ice and the onshore-offshore gradient in pre-winter zooplankton assemblages in southeastern Beaufort Sea.J. Mar. Syst.74994–1011. 10.1016/j.jmarsys.2007.09.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  DicksonA. G.SabineC. L.ChristianJ. R. (2007). Guide to Best Practices for Ocean CO2 Measurements.Sidney BC: North Pacific Marine Science Organization.

  • Google Scholar
  • View reference in article

• 32

  DoubledayA. J.HopcroftR. R. (2015). Interannual patterns during spring and late summer of larvaceans and pteropods in the coastal Gulf of Alaska, and their relationship to pink salmon survival.J. Plankton Res.37134–150. 10.1093/plankt/fbu092

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  EisnerL. B.NappJ. M.MierK. L.PinchukA. I.AndrewsA. G.III (2014). Climate-mediated changes in zooplankton community structure for the eastern Bering Sea.Deep Sea Res. Part II Top. Stud. Oceanogr.109157–171.

  • Google Scholar
  • View reference in article

• 34

  EisnerL. B.PinchukA. I.KimmelD. G.MierK. L.HarpoldC. E.SiddonE. C. (2017). Seasonal, interannual, and spatial patterns of community composition over the eastern Bering Sea shelf in cold years. Part I: zooplankton.ICES J. Mar. Sci.7572–86. 10.1093/icesjms/fsx156

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  EdgarR. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput.Nucleic Acids Res.321792–1797. 10.1093/nar/gkh340

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  EvansW.MathisJ. T.WinsorP.StatscewichH.WhitledgeT. E. (2013). A regression modeling approach for studying carbonate system variability in the northern Gulf of Alaska.J. Geophys. Res. Oceans118476–489. 10.1029/2012JC008246

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  FabryV. J. (1989). Aragonite production by pteropod molluscs in the subarctic Pacific.Deep Sea Res. Part A Oceanogr. Res. Pap.361735–1751. 10.1016/0198-0149(89)90069-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  FeelyR. A.AlinS. R.CarterB.BednaršekN.HalesB.ChanF.et al (2016). Chemical and biological impacts of ocean acidification along the west coast of North America.Estuar. Coast. Shelf Sci.183260–270. 10.1016/j.ecss.2016.08.043

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  FolmerO.BlackM.HoehW.LutzR.VrijenhoekR. (1994). DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates.Mol. Mar. Biol. Biotechnol.3294–299.

  • Google Scholar
  • View reference in article

• 40

  GanneforsC.BöerM.KattnerG.GraeveM.EianeK.GulliksenB.et al (2005). The Arctic sea butterfly Limacina helicina: lipids and life strategy.Mar. Biol.147169–177. 10.1007/s00227-004-1544-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  GarciaH. E.LocarniniR. A.BoyerT. P.AntonovJ. I.BaranovaO. K.ZwengM. M.et al (2013). in World Ocean Atlas 2013, Volume 4: Dissolved Inorganic Nutrients (Phosphate, Nitrate, Silicate), Vol. 76edsLevitusS.MishonovA. (Washington DC: NOAA Atlas NESDIS), 25.

  • Google Scholar
  • View reference in article

• 42

  GascaR.JanssenA. W. (2014). Taxonomic review, molecular data and key to the species of Creseidae from the Atlantic Ocean.J. Molluscan Stud.8035–42. 10.1093/mollus/eyt038

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  GellerJ.MeyerC.ParkerM.HawkH. (2013). Redesign of PCR primers for mitochondrial cytochrome c oxidase subunit I for marine invertebrates and application in all-taxa biotic surveys.Mol. Ecol. Resour.13851–861. 10.1111/1755-0998.12138

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  GibsonG. A.SpitzY. H. (2011). Impacts of biological parameterization, initial conditions, and environmental forcing on parameter sensitivity and uncertainty in a marine ecosystem model for the Bering Sea.J. Mar. Syst.88214–231. 10.1016/j.jmarsys.2011.04.008

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  HaradssonC.AndersonL. G.HassellövM.HulthS.OlssonK. (1997). Rapid high-precision potentiometric titration of alkalinity in ocean and sediment pore waters.Deep Sea Res. Part I442031–2044. 10.1016/s0967-0637(97)00088-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  HauriC.PagèsR.McDonnellA. M. P.StueckerM. F.DanielsonS. L.HedstromK.et al (2021). Modulation of ocean acidification by decadal climate variability in the Gulf of Alaska.Nat. Commun. Earth Environ.(in press). 10.1038/s43247-021-00254-

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  HauriC.SchultzC.HedstromK.DanielsonS.IrvingB.DoneyS. C.et al (2020). A regional hindcast model simulating ecosystem dynamics, inorganic carbon chemistry, and ocean acidification in the Gulf of Alaska.Biogeosciences173837–3857. 10.5194/bg-17-3837-2020

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  HebertP. D.RatnasinghamS.De WaardJ. R. (2003). Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species.Proc. R. Soc. Lond.270(Suppl. 1)S96–S99.

  • Google Scholar
  • View reference in article

• 49

  IPCC (2014). in Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edsCore Writing TeamPachauriR. K.MeyerL. A. (Geneva: IPCC), 151.

  • Google Scholar
  • View reference in article

• 50

  IPCC (2019). in IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edsPörtnerH.-O.RobertsD. C.Masson-DelmotteV.ZhaiP.TignorM.PoloczanskaE.et al (Geneva: IPCC).

  • Google Scholar
  • View reference in article

• 51

  JanssenA. W.BurridgeA. K.MekkesL.PeijnenburgK. T. C. A. (2018). History of the Clio pyramidata-complex (Mollusca, Gastropoda, Euthecosomata, Cliidae) in the northern and central Atlantic Ocean, Caribbean and Mediterranean, restoring its original concept.Basteria8283–102.

  • Google Scholar
  • View reference in article

• 52

  JanssenA. W.BushS. L.BednaršekN. (2019). The shelled pteropods of the northeast Pacific Ocean (Mollusca: Heterobranchia, Pteropoda).Zoosymposia13305–346. 10.11646/zoosymposia.13.1.22

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  JohnsonK. M.AndrewG. D.EischeidG.CatherineG.GuentherP.KeyR. M.et al (1998). Coulometric total carbon dioxide analysis for marine studies: assessment of the quality of total inorganic carbon measurements made during the US Indian Ocean CO2 Survey 1994–1996.Mar. Chem.6321–37.

  • Google Scholar
  • View reference in article

• 54

  KearneyK.HermannA.ChengW.OrtizI.AydinK. (2020). A coupled pelagic–benthic–sympagic biogeochemical model for the Bering Sea: documentation and validation of the BESTNPZ model (v2019.08.23) within a high-resolution regional ocean model.Geosci. Model Dev.13597–650. 10.5194/gmd-13-597-2020

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  KobayashiH. A. (1974). Growth cycle and related vertical distribution of the thecosomatous pteropod Spiratella (“Limacina”) helicina in the central Arctic Ocean.Mar. Biol.26295–301. 10.1007/BF00391513

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  KumarS.StecherG.TamuraK. (2016). MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets.Mol. Biol. Evol.331870–1874. 10.1093/molbev/msw054

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  LalliC. M.GilmerR. W. (1989). Pelagic Snails: The Biology of Holoplanktonic Gastropod Mollusks.Palo Alto CA: Stanford University Press.

  • Google Scholar

• 58

  LauvsetS. K.KeyR. M.OlsenA.Van HeuvenS.VeloA.LinX.et al (2016). A new global interior ocean mapped climatology: The 1× 1 GLODAP version 2.Earth Syst. Sci. Data8325–340.

  • Google Scholar

• 59

  LeeK.FeelyR. A. (2021). News and views: dissolution resolution.Nat. Geosci.14356–358. 10.1038/s41561-021-00765-6

  • CrossRef
  • Google Scholar

• 60

  LeeK.KimT. W.ByrneR. H.MilleroF. J.FeelyR. A.LiuY. M. (2010). The universal ratio of boron to chlorinity for the North Pacific and North Atlantic oceans.Geochim. Cosmochim. Acta741801–1811.

  • Google Scholar

• 61

  LeighJ. W.BryantD. (2015). popart: full-feature software for haplotype network construction.Methods Ecol Evol61110–1116. 10.1111/2041-210X.12410

  • CrossRef
  • Google Scholar

• 62

  LewisE.WallaceD. W. R. (1998). Program Developed for CO2 System Calculations, ORNL/CDIAC-105, Carbon Dioxide Inf. Anal. Cent., Oak Ridge Natl. Lab., Oak Ridge, Tenn. 38. https://salish-sea.pnnl.gov/media/ORNL-CDIAC-105.pdf(accessed January 2019).

  • Google Scholar

• 63

  LischkaS.BüdenbenderJ.BoxhammerT.RiebesellU. (2011). Impact of ocean acidification and elevated temperatures on early juveniles of the polar shelled pteropod Limacina helicina: mortality, shell degradation, and shell growth.Biogeosciences (BG)8919–932. 10.5194/bg-8-919-2011

  • CrossRef
  • Google Scholar

• 64

  LongW. C.SwineyK. M.HarrisC.PageH. N.FoyR. J. (2013). Effects of Ocean Acidification on Juvenile Red King Crab (Paralithodes camtschaticus) and Tanner Crab (Chionoecetes bairdi) Growth, Condition, Calcification, and Survival.PLoS One8:e60959. 10.1371/journal.pone.0060959

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 65

  LongW. C.SwineyK. M.FoyR. J. (2016). Effects of high pCO2 on Tanner crab reproduction and early life history, Part II: carryover effects on larvae from oogenesis and embryogenesis are stronger than direct effects.ICES J. Mar. Sci.73836–848. 10.1093/icesjms/fsv251

  • CrossRef
  • Google Scholar

• 66

  LuekerT. J.DicksonA. G.KeelingC. D. (2000). Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K1 and K2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium.Mar. Chem.70105–119. 10.1016/s0304-4203(00)00022-0

  • CrossRef
  • Google Scholar

• 67

  MaasA. E.Blanco-BercialL.LawsonG. L. (2013). Reexamination of the species assignment of diacavolinia pteropods Using DNA barcoding.PLoS One8:e53889. 10.1371/journal.pone.0053889

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  MackasD. L.GalbraithM. D. (2012). Pteropod time-series from the NE Pacific.ICES J. Mar. Sci.69448–459. 10.1093/icesjms/fsr163

  • CrossRef
  • Google Scholar

• 69

  MajewskiA. R.WalkuszW.LynnB. R.AtchisonS.EertJ. (2016). Distribution and diet of demersal Arctic Cod, Boreogadus saida, in relation to habitat characteristics in the Canadian Beaufort Sea.Polar Biol.391087–1098. 10.1007/s00300-015-1857-y

  • CrossRef
  • Google Scholar

• 70

  MannoC.PeckV. L.TarlingG. A. (2016). Pteropod eggs released at high pCO 2 lack resilience to ocean acidification.Sci. Rep.6:25752. 10.1038/srep25752

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  MathisJ. T.CrossJ. N.BatesN. R. (2011). Coupling primary production and terrestrial runoff to ocean acidification and carbonate mineral suppression in the eastern Bering Sea.J. Geophys. Res. Oceans11610.1029/2010JC006453

  • CrossRef
  • Google Scholar

• 72

  MathisJ. T.CrossJ. N.MonacciN.FeelyR. A.StabenoP. (2014). Evidence of prolonged aragonite undersaturations in the bottom waters of the southern Bering Sea shelf from autonomous sensors.Deep Sea Res. Part II Top. Stud. Oceanogr.109125–133. 10.1016/j.dsr2.2013.07.019

  • CrossRef
  • Google Scholar

• 73

  MathisJ. T.CrossJ. N.EvansW.DoneyS. C. (2015). Ocean acidification in the surface waters of the Pacific-Arctic boundary regions.Oceanography28122–135. 10.5670/oceanog.2015.36

  • CrossRef
  • Google Scholar

• 74

  McGowanJ. A. (1963). “Geographical variation in Limacina helicina in the North Pacific,” in Speciation in the Sea, edsHardingJ. P.TebbleN. (London: Publications of the Systematics Association), 109–128.

  • Google Scholar

• 75

  MekkesL.RenemaW.BednaršekN.AlinS. R.FeelyR. A.HuismanJ.et al (2021). Pteropods make thinner shells in the upwelling region of the California current ecosystem.Sci. Rep.11:1731.

  • Google Scholar

• 76

  NiemiA.BednaršekN.MichelC.FeelyR. A.WilliamsW.Azetsu-ScottK. (2021). Widespread severe pteropod shell dissolution in Amundsen Gulf.Front. Mar. Sci.8:600184. 10.3389/fmars.2021.600184

  • CrossRef
  • Google Scholar

• 77

  OakesR. L.SessaJ. A. (2020). Determining how biotic and abiotic variables affect the shell condition and parameters of Heliconoides inflatus pteropods from a sediment trap in the Cariaco Basin.Biogeosciences171975–1990. 10.5194/bg-17-1975-2020

  • CrossRef
  • Google Scholar

• 78

  PeckV. L.TarlingG. A.MannoC.HarperE. M.TynanE. (2016). Outer organic layer and internal repair mechanism protects pteropod Limacina helicina from ocean acidification.Deep Sea Res. Part II Top. Stud. Oceanogr.12741–52. 10.1016/j.dsr2.2015.12.005

  • CrossRef
  • Google Scholar

• 79

  PeckV. L.OakesR. L.HarperE. M.MannoC.TarlingG. A. (2018). Pteropods counter mechanical damage and dissolution through extensive shell repair.Nat. Commun.9:264.

  • Google Scholar

• 80

  PetersW.JacobsonA. R.SweeneyC.AndrewsA. E.ConwayT. J.MasarieK.et al (2007). An atmospheric perspective on North American carbon dioxide exchange: carbontracker.Proc. Natl. Acad. Sci.U.S.A.10418925–18930. 10.1073/pnas.0708986104

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 81

  PilcherD. J.NaimanD. M.CrossJ. N.HermannA. J.SiedleckiS. A.GibsonG. A.et al (2019). Modeled effect of coastal biogeochemical processes, climate variability, and ocean acidification on aragonite saturation state in the Bering Sea.Front. Mar. Sci.5:508. 10.3389/fmars.2018.00508

  • CrossRef
  • Google Scholar

• 82

  SahaS.MoorthiS.PanH. L.WuX.WangJ.NadigaS.et al (2010). The NCEP climate forecast system reanalysis.Bull. Am. Meteorol. Soc.911015–1057. 10.1175/2010BAMS3001.1

  • CrossRef
  • Google Scholar

• 83

  SahaS.MoorthiS.WuX.WangJ.NadigaS.TrippP.et al (2014). The NCEP climate forecast system version 2.J. Clim.272185–2208. 10.1175/JCLI-D-12-00823.1

  • CrossRef
  • Google Scholar

• 84

  ShchepetkinA. F.McWilliamsJ. C. (2005). The regional oceanic modeling system (ROMS): a split explicit, free-surface, topography-following-coordinate oceanic model.Ocean Modell.9347–404. 10.1016/j.ocemod.2004.08.002

  • CrossRef
  • Google Scholar

• 85

  ShimizuK.KimotoK.NoshitaK.WakitaM.FujikiT.SasakiT. (2018). Phylogeography of the pelagic snail Limacina helicina (Gastropoda: Thecosomata) in the subarctic western North Pacific.J. Molluscan Stud.8430–37. 10.1093/mollus/eyx040

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 86

  ShimizuK.NoshitaK.KimotoK.SasakiT. (2021). Phylogeography and shell morphology of the pelagic snail Limacina helicina in the Okhotsk Sea and western North Pacific.Mar. Biodivers.51:22.

  • Google Scholar

• 87

  SromekL.LasotaR.SzymelfenigM.WolowiczM. (2015a). Genetic evidence for the existence of two species of the “Bipolar” Pelagic Mollusk Clione limacina.Am. Malacol. Bull.33118–120. 10.4003/006.033.0108

  • CrossRef
  • Google Scholar

• 88

  SromekL.LasotaR.WolowiczM. (2015b). Impact of glaciations on genetic diversity of pelagic mollusks: antarctic Limacina antarctica and Arctic Limacina helicina.Mar. Ecol. Prog. Ser.525143–152. 10.3354/meps11237

  • CrossRef
  • Google Scholar

• 89

  StackpooleS.ButmanD.ClowD.VerdinK.GagliotiB.StrieglR. (2016). “Chapter 8. carbon burial, transport, and emission from inland aquatic ecosystems in Alaska,” in Baseline and Projected Future Carbon Storage and Greenhouse-Gas Fluxes in Ecosystems of Alaska, edsZhuZ.DavidA.GuireMc (Reston VA: U.S. Geological Survey).

  • Google Scholar

• 90

  StackpooleS. M.ButmanD. E.ClowD. W.VerdinK. L.GagliotiB. V.GenetH.et al (2017). Inland waters and their role in the carbon cycle of Alaska.Ecol. Appl.271403–1420. 10.1002/eap.1552

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 91

  StockC. A.DunneJ. P.JohnJ. G. (2014). Global-scale carbon and energy flows through the marine planktonic food web: an analysis with a coupled physical–biological model.Prog. Oceanogr.1201–28. 10.1016/j.pocean.2013.07.001

  • CrossRef
  • Google Scholar

• 92

  SulpisO.JeanssonE.DinauerA.LauvsetS. K.MiddelburgJ. J. (2021). Calcium carbonate dissolution patterns in the ocean.Nat. Geosci.14423–428. 10.1038/s41561-021-00743-y

  • CrossRef
  • Google Scholar

• 93

  TsujinoH.UrakawaS.NakanoH.SmallR. J.KimW. M.YeagerS. G.et al (2018). JRA-55 based surface dataset for driving ocean–sea-ice models (JRA55-do).Ocean Modell.13079–139. 10.1016/j.ocemod.2018.07.002

  • CrossRef
  • Google Scholar

• 94

  TsurumiM.MackasD. L.WhitneyF. A.DiBaccoC.GalbraithM. D.WongC. S. (2005). Pteropods, eddies, carbon flux, and climate variability in the Alaska gyre.Deep Sea Res. Part II Top. Stud. Oceanogr.521037–1053. 10.1016/j.dsr2.2005.02.005

  • CrossRef
  • Google Scholar

• 95

  van OostendeN.DussinR.StockC. A.BartonA. D.CurchitserE.DunneJ. P.et al (2018). Simulating the ocean’s chlorophyll dynamic range from coastal upwelling to oligotrophy.Prog. Oceanogr.168232–247. 10.1016/j.pocean.2018.10.009

  • CrossRef
  • Google Scholar

• 96

  WalkuszW.WilliamsW. J.KwasniewskiS. (2013). Vertical distribution of mesozooplankton in the coastal Canadian Beaufort Sea in summer.J. Mar. Syst.12726–35. 10.1016/j.jmarsys.2012.01.001

  • CrossRef
  • Google Scholar

• 97

  WilliamsC. R.DittmanA. H.McElhanyP.BuschD. S.MaherM. T.BammlerT. K.et al (2019). Elevated CO2 impairs olfactory-mediated neural and behavioral responses and gene expression in ocean-phase coho salmon (Oncorhynchus kisutch).Global Chang. Biol.25963–977. 10.1111/gcb.14532

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 98

  WilliamsS. E.ShooL. P.IsaacJ. L.HoffmannA. A.LanghamG. (2008). Towards an integrated framework for assessing the vulnerability of species to climate change.PLoS Biol.6:2621–2626. 10.1371/journal.pbio.0060325

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 99

  ZavolokinA. V.EfimkinA. Y.SlabinskiyA. M.KosenokN. S. (2007). Food supply and trophic relationships of Pacific salmon (Oncorhynchus spp.) and Atka mackerel (Pleurogrammus monopterygius) in the western Bering Sea in fall 2002–2004.N. Pac. Anadr. Fish. Comm. Bull.4127–131.

  • Google Scholar

• 100

  ZavolokinA. V.KulikV. V.ZavarinaL. O. (2014). The food supply of the Pacific salmon of the genus Oncorhynchus in the northwestern Pacific Ocean. 2. comparative characterization and general state.Russ. J. Mar. Biol.40199–207. 10.1134/S1063074014030092

  • CrossRef
  • Google Scholar