Greater exposure of nearshore habitats in the Bering Sea makes fish early life stages vulnerable to climate change

Marine fish species with complex life cycles undergo ontogenetic changes in their physiological and habitat requirements. Therefore, predicting how warming ocean conditions are likely to impact fish populations requires size-(age-) and habitat-specific analyses. We determined the habitat exposure and scope of adaptability of 4 shelf-oriented flatfish species in the Bering Sea to projected climate-driven warming. We quantified present day and end-of-century habitat ex posure based on hindcasts and forecasts of regional ocean circulation models and quantified fish adaptability based on changes in depth distribution and from published thermal tolerances of northern rock sole Lepidopsetta polyxystra, yellowfin sole Limanda aspera, Alaska plaice Pleuro nectes quadrituberculatus, and flathead sole Hippoglossoides elassodon. These 4 species complete their life cycle within the Bering Sea shelf but have different depth preferences and thermal tolerances throughout ontogeny. We found that species or size ranges that occupy the inner shelf, such as northern rock sole, yellowfin sole, and Alaska plaice, are exposed to higher seasonal variability compared to outer shelf species. While these inner shelf species are likely adapted to large seasonal changes in temperature, the future range of seasonal variability was projected to exceed their thermal tolerances. Therefore, we expect species that reside inshore during part of their life cycle and have high temperature sensitivity and limited mobility to be particularly vulnerable to climate change.


INTRODUCTION
Marine fish species with complex life cycles undergo ontogenetic changes in their physiological and habitat requirements (Bartolino et al. 2011). Typically, settled individuals of many shelf or slopeoriented species move from shallow to deep habitats as they age, and coastal and shallow regions are im -portant nursery habitats during juvenile stages (Gibson et al. 2002, Ryer et al. 2010). Therefore, modeling species−environment relationships and predicting how warming ocean conditions are likely to impact fish populations requires size-(age-) (e.g. Barbeaux & Hollowed 2018) and habitat-specific analyses.
During early life stages, fish are limited in their ability to move away from unfavorable conditions, making them more vulnerable to extreme environmental changes (Llopiz et al. 2014, Dahlke et al. 2020). In addition, the exposure of different fish habitats to environmental change varies greatly within the same ecosystem. While nearshore fish species may be adapted to large seasonal shifts in water temperature, the projected variability may exceed the thermal window that these species can tolerate, particularly for stenothermic species of Arctic and sub-Arctic origins (Mueter & Litzow 2008).
In this study, we determined the vulnerability of 4 shelf-oriented flatfish species in the Bering Sea (BS). Species vulnerability to climate change is influenced by both exposure to different levels of (or variations in) environmental parameters (e.g. temperature) and sensitivity to such exposure (IPCC 2014). Sensitivity is, in turn, a function of adaptive reactions, which include behavioral (such as changes in distribution), physiological (such as acclimation), and microevolutionary responses (Donelson et al. 2019). Here, we focused the sensitivity assessment of behavioral responses by determining the degree to which species change their depth distribution in relation to changes in temperature. We quantified current and future (end of century) habitat exposure based on hindcast and forecast of ocean circulation models tested for the region.
We focused our study on flatfish species that complete their life cycle on the shelf (< 200 m depth) but have varying degrees of distribution and exhibit habitat changes throughout ontogeny: northern rock sole Lepidopsetta polyxystra, which completes its life cycle in the eastern and warmer portion of the BS inner shelf; yellowfin sole Limanda aspera, which has pronounced ontogenetic and seasonal migrations between the western and coldest position of the inner and outer shelf (Wilderbuer et al. 1992); Alaska plaice Pleuronectes quadrituberculatus, which completes its life cycle in the inner and middle central portion of the shelf; and flathead sole Hippoglossoides elassodon, which is an outer shelf species.
In the BS, seasonal and interannual variability in the temperatures of the bottom waters that all 4 of these flatfish species inhabit are tightly coupled to variability in sea ice extent and timing (Stabeno et al. 2001). In the future, bottom temperatures may experience changes in variability due to the effects of climate change on sea ice dynamics. The 4 focal species are relatively less studied compared to other larger and commercially important flatfish species (e.g. Pacific halibut, arrowtooth flounder) and have restricted depth preferences and specialized diets which may limit the ability of individuals to shift their depth dis-tribution (Swartzman et al. 1992, McConnaughey & Smith 2000. Additionally, these species are known to use different habitats from early to adult stages and are, therefore, likely to experience different amounts of variation in environmental conditions during their life cycle (Wilderbuer et al. 1992, Porter & Ciannelli 2018, Nichol et al. 2019. While our analyses focused on the BS, the assessment of species climate vulnerability is germane to many coastal marine systems that are experiencing progressive warming (Holsman et al. 2017).

MATERIALS AND METHODS
The study area included the eastern BS shelf, delimited by the date line to the west and the Alaska Peninsula to the east, and the 30 and 200 m isobaths to the north and south (see Fig. 1). During summer, the eastern BS shelf can be separated into 3 oceanographic domains: (1) the inner shelf, with depth < 50 m and mixed water column; (2) the middle shelf, between 50−100 m depth, characterized by 2 surface fronts, cold bottom temperature, and progressive water column stratification; and (3) the outer shelf, between 100−200 m depth, with stratified water column and influx of bottom warm water from the slope and BS basin. During summer, much of the middle shelf is characterized by a frigid water mass, the cold pool, formed as a result of the previous winter sea ice formation and onset of summer water stratification (Stabeno et al. 2012).
We used catch per unit effort (CPUE), depth, and bottom water temperature data from the National Marine Fisheries Service BS summer bottom trawl surveys from 2000−2018 (Stevenson & Lauth 2019). We chose to focus the analyses on these 2 decades because they include stanzas of warm and cold years and encompass the most extreme variation in bottom water temperatures on the BS shelf (Stabeno et al. 2017). Further, this period includes the most reliable and extensive collection of fish length data.
Survey start dates ranged from 24 May−12 June over the history of the survey (1982−present), and end dates ranged from 12 July−11 September; the mean (across stations) standard deviation (across years) in sampling date at any specific station in the primary survey polygon was 6.61 d over this period. Bottom temperatures in the middle and outer shelf regions were very stable across these time periods due to strong stratification. Bottom temperatures in the inner domain can vary more widely across these time scales due to surface heating, but interannual variability still tended to be higher than within-year variability during the sampling period (see Kearney 2021 for details). For most years, the spatial ex tent of the survey spanned from the Alaska Peninsula to approximately 61°N, and the depths at trawl stations ranged from 30−192 m (see Fig. 1). In 2010 and 2017, the survey expanded into the northern BS shelf; however, we excluded this area from our analysis.
We used measurements from the survey to quantify average temperatures and coefficients of variation (CVs) ( Table 1). However, given the 2 mo duration of the survey, these measurements are not synoptic and may compound spatial and temporal variability of the bottom temperature on the shelf. Therefore, regional model output was used in place of the survey temperatures to quantify habitat exposure. Survey temperatures were used to classify years as warm or cold based on whether the summer average temperature of the middle shelf (50−100 m depth) was above or below the reference average from 2000−2018 (Table 1). Identification of extreme marine heatwave years was based on Carvalho et al. (2021).
We used a high-resolution regional ocean model to quantify bottom temperature magnitude, range, and variability (i.e. habitat exposure) for the southeastern shelf region across both seasonal and interannual periods. The model, known as Bering10K, is an implementation of the Regional Ocean Modeling System (ROMS), a free-surface, primitive equation hydrographic model (Shchepetkin & McWilliams 2005, Haid vogel et al. 2008. Its domain spans the BS and northern Gulf of Alaska with 10 km horizontal resolution and 30 terrain-following depth levels. Bering10K has demonstrated high skill in capturing spatiotemporal variability in bottom temperature across the BS region (Kearney et al. 2020, Kearney 2021. We applied our exposure analysis to 2 sets of Bering10K simulations. The first was a hindcast simulation driven by reanalysis-based atmospheric and lateral ocean boundary conditions (see Kearney et al. 2020 for full details). The full simulation spans 1970− 2021; we limited our analysis to the period of 1990− 2020, which encompasses typical decadal temperature variability and overlaps with the groundfish survey data. We also looked at a long-term multimodel forecast driven by global-scale models from the sixth phase of the Climate Model Intercomparison Project (CMIP6). Three specific parent models (MIROC ES2L, CESM2, and GFDL ESM4) were chosen because they encompass the envelope of variability in long-term trends within the larger CMIP6 suite. The CMIP6 global earth system models themselves do not fully capture biophysically relevant features, such as the formation of < 2°C bottom water on the shelf, due to their coarse horizontal resolution (~1°) (see Kearney et al. 2020 for more details). Dynamically downscaling (i.e. driving a highresolution regional ocean model using atmospheric and lateral ocean boundary condition data from the CMIP6 models) allows us to simulate important regional-scale processes with higher resolution while preserving the larger-scale climate state predicted by the global models. Unfortunately, we were unable to downscale the entire CMIP6 suite due to computational and data storage limitations, hence the choice of specific ensemble members from within the larger collection to represent the uncertainty envelope.
We quantified variability for 3 periods within these forecasts: current (1985−2015), mid-century (2035− 2065), and end-of-century (2070−2100). While it is somewhat redundant to show results for both the hindcast simulation and the historical periods of the forecasts, we chose to do so because the hindcast simulation is driven by a data-assimilating reanalysis and therefore can be used for direct interannual comparison. Also, the global climate model data were not bias-corrected or otherwise adjusted to observations for these particular downscaling simulations. Therefore, the 1990−2020 portion of each climate model projection includes any biases or other mismatches that are present in these parent models. to BS bottom temperature, these mismatches are small but not negligible. When comparing the downscaled earth system model-forced historical simulations to the reanalysis-forced hindcast, the GFDL and MIROC bottom temperatures showed a small cool bias across the eastern BS shelf region while the CESM bottom temperatures were biased a bit warm. The mid-century and end-of-century changes were therefore assessed relative to their respective historical time period values rather than using the hindcast values as a baseline. Variability was measured via the following yearly bottom temperature metrics, calculated for each individual 10 km grid cell for months between May and October (when juvenile flatfish stages are settled in their nursery areas): (1) monthly mean temperature: mean temperature for dates falling within each month; (2) maximum weekly mean temperature: max. temperature between 1 May and 31 Oct; (3) minimum weekly mean temperature: min. temperature between 1 May and 31 Oct; (4) temperature range (maximum − minimum): difference between 2 and 3.
For each metric, we calculated the climatological mean across each 30 yr period and the CV, defined as the standard deviation of the metric across the 30 yr period divided by the climatological mean (in K). Calculations were applied to weekly averaged output. The forecast metrics were then averaged across the 3 parent models to acquire a single multi-model forecast value. The 8 metrics (4 means and corresponding CVs) were used to define magnitude, range, and respective variation of bottom temperature at each ROMS grid location, and allowed us to differentiate between interannual and seasonal temperature variability. Specifically, the temperature range de fines locations with the greatest seasonal temperature va riability, while the corresponding interannual CV defines locations with the greatest interannual variability.
To quantify species adaptability, we assessed changes in distribution and temperature exposure using length and abundance data from the trawl survey. Length data represent a subsample of the haul. The proportion of individuals of a given sex and length in the subsample was multiplied by the number of fish of that species in the haul to provide an estimate of the total number of fish by length and sex in each haul. For each species, we quantified the weighted (by CPUE) mean depth and temperature (from survey) across 8 length quantiles. The estimated means and variability for each length quantile are shown as boxplots to compare cold:warm years and cold years:extreme marine heatwave years (Table 1). Size-specific image plots of individual density in relation to either bottom temperature or depth of occupancy were calculated using 2-dimensional kernel density smoothing and are plotted against length values for each species, comparing warm and cold groupings using scaled-up length as a continuous variable (i.e. no binning).

Metrics of habitat exposure using the Bering10K model
During the period 1990−2020, the simulated mean summer (July) bottom temperature of the BS shelf reproduced the known signature for the area, with a cold water mass in the middle shelf ranging from −1.8 to 3.0°C enclosed between warmer water masses in the inner and outer shelf ranging from 3.0−12.0 and 3.0−6.0°C, respectively (Fig. 1). The spatial patterns of mean bottom temperature and variation from survey data corroborated the patterns in the ROMS model hindcast simulations during the summer months (see Fig. S1 in the Supplement at www.int-res.com/ articles/ suppl/ m684 p091_supp.pdf). This spatial temperature signature was replicated through May− October, albeit with variation in the temperature values (Fig. S2).
The interannual CVs of mean July temperatures were highest in the shallow inner and north portions of the shelf and lowest along the deep outer shelf and the northwesternmost portion of the inner and middle shelf (Fig. 1). In August and September, peaks of temperature variation occur offshore, toward the middle shelf. The seasonal May through October average range of bottom temperatures is highest in the northeast inner shelf, where it reaches values of 10−11°C and where the maximum seasonal bottom temperature can reach 13°C. In contrast, the May− October temperature range was less than 3°C in the outer shelf region, denoting a lower exposure of these habitats to seasonal temperature variations (Fig. S2). However, the interannual CV of temperature range was higher in the outer shelf than in the inner shelf, indicating that despite lower seasonal variability, deeper habitats are exposed to high interannual variability of bottom temperature (Fig. S2).
The medium (2035−2065) and long-term (2070− 2100) forecasts of bottom temperature predict a warmer BS shelf due to an increase of both May− Oct minimum and maximum temperatures throughout the 3 domains (inner, middle, and outer shelf; Figs. S3 & S4). During the historical record (1990− 2020), the model predicts an average of 0.32 wk yr −1 when the bottom temperature of the inner shelf reaches values of 12−13°C. For the inner shelf, the forecast models predict an increase of between 0.36 (GFDL) and 5.96 (CESM) wk yr −1 of temperatures over 13°C by the end of the century relative to the current reference period. Temperatures in the middle and outer shelf are also expected to increase by 3−4°C, but their maximum values will continue to be lower compared to the inner shelf (Figs. S3 & S4).

Ontogenetic patterns of species distributions and exposure
The distribution of rock sole (Fig. 2) is concentrated in the northeast reach of the surveyed area. As they age and grow, they occur in deeper areas of the eastern shelf, with the exception of the largest individuals concentrated in shallower areas. An isolated and stable cluster of abundance is localized around the Pribilof Islands, in the outer shelf. Rock sole average depth distribution was around 40 m for the youngest and smallest individuals captured in the survey (40− 150 mm) and gradually increased to 50−60 m in larger individuals, with the largest individuals examined moving back to shallower waters. These depth distribution patterns did not change between cold and warm years. In these habitats, the fish are exposed to temperature variations averaging 2.5°C, especially during early life stages (Fig. 2).
Yellowfin sole (Fig. 3) is also an inner-middle shelf species. During early and smaller life stages, these fish are almost exclusively found in the cold and northwestern portion of the sampled region. As they grow, they occur in the deeper areas of the middle shelf and the variability around the mean depth distribution increases, as underscored by the expansion of spatial occupancy of intermediate and larger individuals. Their average depth distribution was concentrated around 30−35 m for the youngest and smallest individuals examined (20−180 mm) and reached a maximum of 45−50 m at intermediate and larger size ranges. Noticeably, the variability around the mean depth distribution increased with size, as highlighted by the expansion of spatial occupancy of intermediate and larger individuals. These depth distribution patterns did not change between cold and warm years, and the average range of temperature exposure between warm and cold years varied by about 1.5−2.0°C (Fig. 3).
Alaska plaice has a distinctive middle shelf distribution throughout most of its life cycle (Fig. 4). Only the smaller and younger life stages are found on the inner shelf. Their average depth distribution ranged  During years with marine heatwaves, the 4 species are exposed to markedly increased temperatures (Figs. S5−S8). For example, the smaller juvenile stages of northern rock sole experience an average change of water temperature of about 4−5°C. Smaller and younger stages of yellowfin sole (Fig. S6) and Alaska plaice (Fig. S7) experience an increase of about 3−4°C. Larger flathead sole (Fig. S8) experience higher water temperatures during heatwave years than during cold years.

DISCUSSION
Species adaptability to climate change is a multipronged process that involves behavioral, physiological, and genetic responses (Visser 2008, Donelson et al. 2019. For marine fish, behavioral responses can occur when individuals move from less suitable toward more suitable areas. The rate at which these responses occur varies not only across species and marine systems (Cheung et al. 2009) but also over life history stages (Ciannelli et al. 2022). During ontogenetic development, marine fish species occupy different habitats and exhibit different vulnerabilities to climate change (Petitgas et al. 2013). For example, flatfish species that complete their life cycle on the shelf typically move from shallow coastal habitats during early life stages toward deeper habitats as individuals age and mature (Gibson et al. 2002). Our study demonstrates that flatfish stages inhabiting nearshore areas are exposed to double jeopardy from climate change. First, these stages are exposed to higher seasonal variability than those occupying deeper water. Second, their limited ability to move deeper when temperatures increase makes them more vulnerable to climate change.
While shallow water species/stages are likely physiologically adapted to large changes in temperature, the limits of their adaptability will be tested when the seasonal envelope of exposure exceeds their thermal optimum. For example, in the BS, the thermal optimum for juvenile rock sole is ~13°C, and survival starts to rapidly decline beyond that threshold (Ryer et al. 2012). The average temperature of the BS inner shelf habitats in 2070−2100 is expected to surpass 13°C for stretches of 1 mo or longer. Other inner shelf species such as yellowfin sole and Alaska plaice that have a lower thermal optimum and higher thermal sensitivity than rock sole (T. Hurst pers. comm.) are likely to be more vulnerable to extreme heat waves and continuous warming over the next few decades.
The BS shelf is physically and biologically heterogeneous, hosting a variety of habitats and species assemblages, mostly organized along bathymetry and temperature gradients (Baker & Hollowed 2014, Sigler et al. 2017. During summer, the bottom temperatures of the middle and outer shelf are less variable than that of the inner shelf. During late summer and fall, bottom temperature becomes more variable offshore. In general, offshore and deeper habitats are exposed to lower seasonal temperature variability during summer but greater interannual variability compared to shallower habitats of the BS. However, the maximum temperatures of the inner shelf are projected to increase beyond the optimal temperatures for growth of stenotherm species of Arctic and sub-Arctic origins, while middle and outer shelf temperatures will continue to remain below these optimal temperatures. Therefore, we expect that species and stages that reside inshore and that cannot move to other habitats are likely to be more vulnerable to climate change. Apart from flathead sole, all species examined in this study spend their early life stage in some portion of the BS inner shelf, and yellowfin sole early stages are almost exclusively found in shallow nearshore habitats. Flathead sole stands apart from all other species: they live in deeper and less thermally exposed habitats and move to shallower regions as they grow. This is potentially a good strategy to cope with increasing water temperatures. A species' sensitivity to an environmental variable is proportional to the slope of the relationship between that variable and a particular trait (reaction norm ;Visser 2008). It follows that the slope of the relationship between temperature and individual growth is a measure of the species' thermal sensitivity -with a more positive correlation indicating a higher sensitivity to temperature changes. Lower sensitivity to changes in temperature could be due to a variety of reasons, including physiological or behavioral responses or less exposure to thermal changes, if temperature is measured as an index over different habitats (e.g. at regional scales). Matta et al. (2010) found that otolith growth increments of Alaska plaice and yellowfin sole were more positively correlated with the summer water temperature of the eastern BS shelf area than that of northern rock sole, underscoring the thermal sensitivity and potential vulnerability to extreme high temperatures of the former 2 species. Northern rock sole was the least correlated, probably because of their ability to cope with higher temperature and their relatively higher habitat occupancy, especially during high abundance years, compared to the other 2 species examined.
Other variables, including the presence of infaunal and epibenthic prey, sediment size and lithology, predation, and fishing pressure are known to affect flatfish distribution. In previous studies conducted in the BS, several of these variables contributed to spatially anchoring species to specific habitats (Mc -Connaughey & Smith 2000, Spencer 2008). In terms of sediment type, inshore-oriented species occur within a mixture of sand and gravel, while offshore species are associated with fine sand and mud (Spencer 2008). Organic content and sediment size are correlated with higher trophic level community composition and abundance in benthic habitats (Greb meier et al. 1989), and bathymetry is negatively correlated with sediment size. Not surprisingly, most of the shelf flatfish species in the BS have a strong affinity with depth (Swartzman et al. 1992). Spencer (2008) examined changes in distribution of several flatfish species, including the 4 studied here, in relation to the BS cold pool extent. They found that among 6 flatfish species examined, distribution shifts in the southeast region of the shelf were only correlated to the extent of the cold pool in flathead sole and northern rock sole. Bartolino et al. (2011) found that yellowfin sole juvenile stages shift their distribution in relation to water temperature, but to a lesser extent than adult individuals. In a follow-up study (Bartolino et al. 2012), the authors examined the effect of commercial harvest and found evidence of local depletion. Results from other systems confirm that fishing, in addition to temperature, is an important driver of fish distribution (Ciannelli et al. 2013, Frank et al. 2018. For example, Engelhard et al. (2011) found that the distribution of sole Solea solea in the North Sea was affected by fishing pressure. In addition, both sole and plaice Pleuronects platessa changed their distribution in relation to climate variability, albeit in opposite directions -sole moved southward and plaice moved northward.
There is mounting evidence of climate-driven shifts in species distributions and vulnerability in marine ecosystems worldwide (Dulvy et al. 2008, Pinsky et al. 2019, Gervais et al. 2021. Moreover, a growing number of studies are documenting the varied responses of marine organisms to climate change impacts owing to differences in physiological tolerances among taxonomic groups, types of performance measured, and life stages (Harvey et al. 2013, Twiname et al. 2019, Dahlke et al. 2020. However, our understanding of the vulnerability of marine species and the implications for their geographic distribution and broader ecosystem-level processes under a changing climate is still very limited. Here, we illustrated that within a large and oceanographically diverse region, such as the eastern BS shelf, exposure varies greatly across habitat types, and species sensitivity changes across ontogeny, underscoring the importance of both habitat-specific and size-or age-specific assessments. Results from this study can inform how we represent exposure and sensitivity in climate vulnerability assessments (e.g. Holsman et al. 2017). In previous assessments conducted for fish and invertebrate stocks in the BS, species were found to have limited exposures to climate change (Spencer et al. 2019); however, the authors only carried the projections to 2039. Our initial assessment also highlights the need for more comprehensive analyses of fish sensitivity across ontogeny, such as lab work to better understand physiological limits and thus how fish life stages will be impacted under future climate scenarios.