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

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.

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 shellbuilding 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

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(IPCC, , 2019Carter et al., 2017Carter et al., , 2019Bindoff 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(Long et al., , 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;Lischka et al., 2011;Bednaršek et al., 2014Bednaršek et al., , 2019Feely 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(Eisner et al., , 2017Janssen 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(Eisner et al., , 2017(Eisner et al., , 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 (sensu Williams 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 , 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., 2018Shimizu et al., , 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.

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).

Pteropod Sampling
The list of stations where pteropods were collected, along with the data on abundance (reported in individuals per m 2 ), 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 2 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)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007), as described in Stock et al. (2014). The model was run from 1980 to 2013 and was forced with atmospheric pCO 2 observations from the Mauna Loa CO 2 record ( 1 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 1 https://www.esrl.noaa.gov/gmd/ccgg/trends/data.html TA concentrations in the freshwater were based on available observations (Stackpoole et al., 2016(Stackpoole et al., , 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 shipbased 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 2 varied seasonally following output from the NOAA Earth System Frontiers in Marine Science | www.frontiersin.org

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.

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 (Xray focus spot diameter, 0.8 µm; X-ray tube voltage, 80 kV; detector array size, 1024 × 1024; 1800 projections/360 • , fourtimes 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: (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 (sensu Bednarš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(Bednaršek et al., , 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.

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 2 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).

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).
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).
Frontiers in Marine Science | www.frontiersin.org

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).

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).

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 2 = 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 2 = 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 2 = 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 2 = 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 2 = 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 2 = 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.

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).
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 morphotypespecific.

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).
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 FIGURE 9 | Transect plots of aragonite saturation state ( ar ; based on the observational data) in the Amundsen Gulf in the Canadian Arctic sampled in August 2018. Transects extend across Dolphin and Union Strait (A); Minto Inlet (B) and across the mouth of the Gulf from Cape Bathurst (C). Red boxes on the transect maps include the biological stations investigated for pteropod shell dissolution (previously published in Niemi et al., 2021). 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.
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  Blue colors indicate shorter, and red color indicate longer time intervals since release. Initial release is indicted by a square. The hydrodynamic output from the model hindcasts is described by Coyle et al. (2019) and Kearney et al. (2020).
over the upper 50-100 m of the water column for the winter months can serve as an approximation of the adult exposure . 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 (sensu Bednarš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 (sensu Comeau 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 (sensu McGowan, 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 (sensu Bednarš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 (sensu Manno 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 OArelated 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.

DATA AVAILABILITY STATEMENT
All sequence data are available in Genbank under accession numbers MZ502688-MZ502863 and in the Supplementary Material. The GoA-COBALT model output is publicly available (https://doi.org/10.24431/rw1k43t, Hauri et al., 2020). Carbonate chemistry data from the hydrographic cruise P16N are archived in the NCEI database: https://www.ncei.noaa.gov/data/oceans/ ncei/ocads/metadata/0163182.html. The carbonate chemistry data from the Amundsen Gulf (2018) are publicly available in the Ocean Carbon Data System (OCADS), NCEI Accession number 0238161 (https://www.ncei.noaa.gov/data/oceans/ncei/ ocads/metadata/0238161.html). The Bering10K model output fields and particle tracking output for the Bering Sea and Gulf of Alaska used in this analysis are available online from the Zenodo library (https://doi.org/10.5281/zenodo.5117852). Model source code for the Bering10K ROMS model is available here (https: //github.com/beringnpz/roms-bering-sea). Pteropod abundance data from the GoA and the Bering regions are archived in the NCEI database, along with the hydrographic P16N data; data from the Amundsen Gulf are currently available upon request (Niemi) and will be submitted into the NCEI after the manuscript acceptance. Pteropod shell dissolution data and shell morphometrics and properties are available upon request (Niemi, Bednaršek; pers. communications), followed by the submission into the NCEI database after the manuscript acceptance.

AUTHOR CONTRIBUTIONS
NB wrote the manuscript with contributions from all coauthors. NB and K-AN conceived the research, analyzed the data, and provided visualizations. CH, DP, and AH provided model outputs. NB directed the zooplankton fieldwork. AN provided biological data from the Amundsen Gulf. RF and CM provided carbonate chemistry data. KK analyzed pteropod shell density and thickness and contributed to the findings about the relevance of environmental data. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank Isadora Jimenez-Hidalgo for her technical expertise and Lorenz Hauser and Ana Ramón-Laca for their advice in genetic analyses. We also thank AFSC team, Phyllis Stabeno and Janet Duffy-Anderson. We are grateful to Colleen Harpold and Nissa Ferm to conduct pteropod sampling in the Bering Sea, as well as Peter Proctor for collecting samples for carbonate chemistry in the Bering Sea, and to Shaun Bell to help assist analyses of physical-chemical parameters from Bering Sea cruises. Carbonate chemistry analyses in Beaufort shelf was done in the lab of Kumiko Azetsu-Scott at BIO, funded by ACCASP. Special thank you goes to Jesse Turner who collected pteropod samples on P16, Brendan Carter to collect carbonate chemistry samples, also to Alison MacDonald for making pteropod sampling on P16 feasible. We also thank S. Punshon and J. Charette for carbonate chemistry sample and data analysis for the Amundsen Gulf. We are grateful Arie Janssen with help on distinguishing different pteropod morphotypes.