TY - JOUR AU - Cooper, Daniel, W AB - Abstract Groundfish species in the Bering Sea are undergoing pronounced changes in spatial distribution and abundance due to warming ocean temperatures. The main drivers of interannual variability in this ecosystem are the alternating warm and cold thermal stanzas. Yellowfin sole (Limanda aspera; YFS) and northern rock sole (Lepidopsetta polyxystra; NRS) are commercially-valuable flatfishes in the Bering Sea and are among the most dominant groundfish species there in numbers and biomass. We examined the variability in the spatial distribution and abundance of juvenile NRS and YFS in relation to the ice and temperature conditions associated with warm-cold thermal shifts from 1982 to 2017. The goal was to assess the implications of the fluctuating thermal environment for Bering Sea flatfish production. We found ice cover and bottom temperature indices in the preceding 1 to 3 years to be the best predictors of NRS juvenile distribution. In contrast, these indices were not significantly correlated with YFS juvenile distribution, which could be an artifact of their relatively low availability to sampling. A warm stanza, as the Bering Sea is currently in, is expected to favor high numbers of NRS juveniles and the northward expansion of their distribution. Introduction Rising sea temperature as a result of increasing anthropogenic radiative forcing (Intergovernmental Panel on Climate Change, 2014) is projected to change marine species distributions and assemblages (Overland et al., 2010; Hunt et al., 2011; Fossheim et al., 2015; Britten et al., 2016). High-latitude marine ecosystems are already seeing pronounced changes in the spatial distribution of some economically important groundfish species (Mueter and Litzow, 2008; Hollowed et al., 2013; Stevenson and Lauth, 2019). Ecosystem shifts in the Bering Sea have been studied since the early 2000s (Benson and Trites, 2002; Grebmeier et al., 2006). The Bering Sea connects the North Pacific and the Arctic Ocean. The Bering Sea within the US Exclusive Economic Zone is approximately divided at latitude 60°N into the eastern (also referred to as southeastern in some literature) Bering Sea (EBS) to the south and the northern Bering Sea (NBS) to the north (Sigler et al., 2011; Figure 1). The EBS shelf is typically divided into three depth-associated biophysical domains: inner (≈0–50 m), middle (≈50–100 m), and outer (≈100–200 m) (Kinder and Schumacher, 1981; Coachman, 1986). The most prominent environmental drivers of ecosystem variability on the shelf are the seasonal ice cover and the “cold pool”—the tongue of <2°C bottom water in the middle shelf formed by the spring ice melt that acts as a barrier to some cold-intolerant demersal species (Wyllie-Echeverria and Wooster, 1998; Stabeno et al., 2012). The EBS alternated between periods of warm and cold temperature anomalies that occurred with high year-to-year variability prior to 1999, but since then the periods have lengthened to multiple years (stanzas) (Stabeno et al., 2012; Duffy-Anderson et al., 2017; Stabeno et al., 2017). Generally, in warm stanzas, there is little or no sea-ice cover over the southern shelf in spring and the southward extent of the cold pool on the middle shelf is curtailed; in cold stanzas, there is extensive sea-ice cover in winter and the cold pool is well-developed in summer. Mean currents on the inner shelf are generally westward in cold stanzas and northward in warm stanzas; they are also stronger in cold stanzas (Stabeno et al., 2012). The latest cold stanza was 2006–2013. By late 2013, the EBS has shifted into the present warm stanza (Stabeno et al., 2012; Wood et al., 2015; Duffy-Anderson et al., 2017). Thermal stanzas and the long-term warming trend in the EBS are shifting species distribution, but their effects on production are unclear (Britten et al., 2016). Figure 1. Open in new tabDownload slide Sampling grid of the bottom-trawl survey on the EBS and NBS shelves (divided by bold line). A sampling station is placed approximately at the centre of a cell (square). At St Matthew Island and the Pribilof Islands there are additional stations (circles) at the corner of cells. The dashed lines delineate the 50-, 100-, and 200-m isobaths. The black polygons (dotted lines) enclose sampling strata 20 and 31. Figure 1. Open in new tabDownload slide Sampling grid of the bottom-trawl survey on the EBS and NBS shelves (divided by bold line). A sampling station is placed approximately at the centre of a cell (square). At St Matthew Island and the Pribilof Islands there are additional stations (circles) at the corner of cells. The dashed lines delineate the 50-, 100-, and 200-m isobaths. The black polygons (dotted lines) enclose sampling strata 20 and 31. The wide and flat EBS shelf supports some of the most productive fisheries in the world (Link, 2016). Benthivorous flatfishes dominate the fish biomass in this ecosystem, of which yellowfin sole (Limanda aspera; YFS) and northern rock sole (Lepidopsetta polyxystra; NRS) are the two most abundant and economically valuable species. The EBS YFS fishery is the largest flatfish fishery in the world. Adult populations of both species are concentrated at depths of ≲100 m. Adult YFS overwinter near the shelf margin and migrate onto the inner shelf around and south of Nunivak Island in April–May to spawn and feed (Wakabayashi, 1989; Wilderbuer et al., 2017). Spawning occurs in May–August in nearshore waters of ≲30-m deep between Bristol Bay and at least as far north as Nunivak Island (Nichol and Acuna, 2001) (Figure 1). Pelagic eggs hatch at ≈2.5-mm standard length (SL) (Matarese et al., 2003). Adult NRS also have separate winter to early spring spawning and summer feeding distributions on the EBS shelf (Wilderbuer and Nichol, 2016). Spawning occurs in December–March at depths of 70–140 m (Stark, 2012; Wilderbuer and Nichol, 2012). Eggs are demersal and hatch at ≈3-mm SL. NRS larvae are present in the plankton from March to June along the Alaska Peninsula with the highest abundance in the vicinity of Unimak Island (Lanksbury et al., 2007). Year-class strength of flatfishes is presumably determined in the early life stages between the pelagic egg and benthic settlement (van der Veer et al., 2015). Temperature in the early life stages is critical to recruitment by affecting egg size, larval duration, size at settlement, and extrinsically, the size of suitable nursery habitat (Rätz and Lloret, 2003; Laurel et al., 2014; Fedewa et al., 2017). Under the Maximum Growth/Optimal Feeding hypothesis, growth is solely determined by the temperature in the nursery, if food is not limited (van der Veer and Witte, 1993). Warmer temperature increases juvenile growth rate and improves their condition, thereby decreasing density-dependent mortality, and a large nursery with favourable temperature, prey availability, and situation relative to larval transport pathways would enhance that advantage (van der Veer et al., 2015). There is no evidence that the inner shelf is food-limited for benthivorous flatfishes (Yeung and Yang, 2014, 2018). Recent beam trawl studies (Cooper et al., 2014; Cooper and Nichol, 2016) suggest that habitat usage and densities by early juvenile NRS vary with thermal stanzas. Age 0–1 NRS were concentrated in the northern part of the EBS inner shelf between Cape Newenham and Nunivak Island during a warm year (2003), but the distribution shifted south to near the Alaska Peninsula during 2 cold years (2008 and 2010) (Cooper et al., 2013, 2014). This spatial variation led to the hypothesis that the northern part is occupied only during warm stanzas when currents favour the transport of larvae northward and onshore, and warmer bottom temperatures conditions are more suitable for growth and survival. There is no comparable information on juvenile YFS. In this study, we examine the variability in the spatial distribution and abundance of juvenile NRS and YFS with temperature fluctuations in the Bering Sea from 1982 to 2017, using data primarily from an annual summer bottom trawl survey conducted by the Alaska Fisheries Science Center (AFSC) of the National Oceanic and Atmospheric Administration (NOAA). Key patterns in the juvenile population are related to ice and temperature conditions associated with thermal stanzas. The spatial patterns and population responses of the two benthic foragers are compared for insights into biological and environmental interactions. Based on the apparent trend in Bering Sea temperature, we speculate on the possible implications for adult production in the near future. Methods Bottom trawl survey The indices of juvenile NRS and YFS abundance and spatial distribution (Table 1) for this study were derived from the annual summer (June–August) bottom trawl surveys of the EBS shelf by the AFSC for the assessment of groundfish and invertebrate stocks from 1982 to 2017 (Armistead and Nichol 1993; Bohle and Bakkala, 1994; Conner and Lauth, 2017). An 83–112 eastern otter trawl was used to sample a target of 356 standard stations between the 20- and 200-m isobaths (Figure 1). Each sample consisted of a 30-min tow at 5.6 km h−1 at a station. Fish and invertebrates were sorted and identified. The catch weights and numbers by taxon were estimated and standardized to area swept. Random samples of up to ∼200 individuals of each commercially important species were lengthed, and the lengthed samples were further randomly subsampled to collect sagittal otoliths for ageing (Armistead and Nichol, 1993). Bottom temperature was recorded at each station using expendable bathythermographs in the early years, and depth and temperature data logger mounted on the headrope of the trawl in the recent years. Annual EBS stock assessment reports were produced with survey catch and estimates of biomass and number by age and sex for each managed species (North Pacific Fishery Management Council, 2017a). Bottom temperature was also reported as a principal indicator of EBS ecosystem status (North Pacific Fishery Management Council, 2017b). Table 1. Indices of (a) the spatial distribution and (b) the abundance of NRS and YFS, and (c) the EBS environment. (a) Spatial: Area of 95% presence (area) Estimated area (ha) of minimum convex polygon enclosing all locations excluding 5% of extreme locations farthest away from the centre of spatial distribution Centre of distribution, longitude (cgx), latitude (cgy) Mean geographic location (decimal degree) Global index of collocation (gic) Spatial overlap of populations. Distance between centre of distribution and average distance between individual fish taken at random and independently from each population (0 = completely distinct, 1 = coincide perfectly) Local index of collocation (lic) Co-occurrence at same locations. Can be ≈0 even if gic high, i.e. same range but not exactly the same sites within range (0 = completely distinct, 1 = completely coincide) (b) Abundance: Catch Annual total survey catch (tons), 1982–2017 Number-at-age (NRSi, YFSi) Annual estimated number-at-age (millions) of fish age i < 10, 1982–2017 (c) Environment: Ice cover index anomaly (ICI) Anomaly (based on 1981–2010 average) of average ice concentration for 1 January–31 May within a 2°×2° box (56–58°N, and 163–165°W), 1979–2017 Ice retreat index (IRI) Number of days after March 15 when the average ice concentration within the box is >10% of the total box area, 1979–2017 Wind stress anomalies, November–April (WS_NA), May–June (WS_MJ) Anomaly (based on 1961–2000 average) of along Alaska Peninsula component of wind stress (N·m−2) at Unimak Pass in winter and spring, 1979–2013 Bottom temperature, stratum 20 (BT20), 31 (BT31) Average temperature (°C) recorded by bottom-trawl in survey strata 20 and 31 (Figure 1), 1982–2017 (a) Spatial: Area of 95% presence (area) Estimated area (ha) of minimum convex polygon enclosing all locations excluding 5% of extreme locations farthest away from the centre of spatial distribution Centre of distribution, longitude (cgx), latitude (cgy) Mean geographic location (decimal degree) Global index of collocation (gic) Spatial overlap of populations. Distance between centre of distribution and average distance between individual fish taken at random and independently from each population (0 = completely distinct, 1 = coincide perfectly) Local index of collocation (lic) Co-occurrence at same locations. Can be ≈0 even if gic high, i.e. same range but not exactly the same sites within range (0 = completely distinct, 1 = completely coincide) (b) Abundance: Catch Annual total survey catch (tons), 1982–2017 Number-at-age (NRSi, YFSi) Annual estimated number-at-age (millions) of fish age i < 10, 1982–2017 (c) Environment: Ice cover index anomaly (ICI) Anomaly (based on 1981–2010 average) of average ice concentration for 1 January–31 May within a 2°×2° box (56–58°N, and 163–165°W), 1979–2017 Ice retreat index (IRI) Number of days after March 15 when the average ice concentration within the box is >10% of the total box area, 1979–2017 Wind stress anomalies, November–April (WS_NA), May–June (WS_MJ) Anomaly (based on 1961–2000 average) of along Alaska Peninsula component of wind stress (N·m−2) at Unimak Pass in winter and spring, 1979–2013 Bottom temperature, stratum 20 (BT20), 31 (BT31) Average temperature (°C) recorded by bottom-trawl in survey strata 20 and 31 (Figure 1), 1982–2017 Open in new tab Table 1. Indices of (a) the spatial distribution and (b) the abundance of NRS and YFS, and (c) the EBS environment. (a) Spatial: Area of 95% presence (area) Estimated area (ha) of minimum convex polygon enclosing all locations excluding 5% of extreme locations farthest away from the centre of spatial distribution Centre of distribution, longitude (cgx), latitude (cgy) Mean geographic location (decimal degree) Global index of collocation (gic) Spatial overlap of populations. Distance between centre of distribution and average distance between individual fish taken at random and independently from each population (0 = completely distinct, 1 = coincide perfectly) Local index of collocation (lic) Co-occurrence at same locations. Can be ≈0 even if gic high, i.e. same range but not exactly the same sites within range (0 = completely distinct, 1 = completely coincide) (b) Abundance: Catch Annual total survey catch (tons), 1982–2017 Number-at-age (NRSi, YFSi) Annual estimated number-at-age (millions) of fish age i < 10, 1982–2017 (c) Environment: Ice cover index anomaly (ICI) Anomaly (based on 1981–2010 average) of average ice concentration for 1 January–31 May within a 2°×2° box (56–58°N, and 163–165°W), 1979–2017 Ice retreat index (IRI) Number of days after March 15 when the average ice concentration within the box is >10% of the total box area, 1979–2017 Wind stress anomalies, November–April (WS_NA), May–June (WS_MJ) Anomaly (based on 1961–2000 average) of along Alaska Peninsula component of wind stress (N·m−2) at Unimak Pass in winter and spring, 1979–2013 Bottom temperature, stratum 20 (BT20), 31 (BT31) Average temperature (°C) recorded by bottom-trawl in survey strata 20 and 31 (Figure 1), 1982–2017 (a) Spatial: Area of 95% presence (area) Estimated area (ha) of minimum convex polygon enclosing all locations excluding 5% of extreme locations farthest away from the centre of spatial distribution Centre of distribution, longitude (cgx), latitude (cgy) Mean geographic location (decimal degree) Global index of collocation (gic) Spatial overlap of populations. Distance between centre of distribution and average distance between individual fish taken at random and independently from each population (0 = completely distinct, 1 = coincide perfectly) Local index of collocation (lic) Co-occurrence at same locations. Can be ≈0 even if gic high, i.e. same range but not exactly the same sites within range (0 = completely distinct, 1 = completely coincide) (b) Abundance: Catch Annual total survey catch (tons), 1982–2017 Number-at-age (NRSi, YFSi) Annual estimated number-at-age (millions) of fish age i < 10, 1982–2017 (c) Environment: Ice cover index anomaly (ICI) Anomaly (based on 1981–2010 average) of average ice concentration for 1 January–31 May within a 2°×2° box (56–58°N, and 163–165°W), 1979–2017 Ice retreat index (IRI) Number of days after March 15 when the average ice concentration within the box is >10% of the total box area, 1979–2017 Wind stress anomalies, November–April (WS_NA), May–June (WS_MJ) Anomaly (based on 1961–2000 average) of along Alaska Peninsula component of wind stress (N·m−2) at Unimak Pass in winter and spring, 1979–2013 Bottom temperature, stratum 20 (BT20), 31 (BT31) Average temperature (°C) recorded by bottom-trawl in survey strata 20 and 31 (Figure 1), 1982–2017 Open in new tab In 2010 and 2017, survey coverage was extended into the NBS shelf, thereby adding another 144 stations (Figure 1; Stevenson and Lauth, 2019). The purpose of the extension was to monitor the impacts of diminishing sea ice on the Bering Sea ecosystem, particularly on species distribution (Hollowed et al., 2007; Stevenson and Lauth, 2019). NBS data from these 2 years were included in this study for a perspective of the juvenile flatfish distributions across the entire US Bering Sea shelf. Juveniles were relatively rare in bottom trawl survey catch until about age 5, or ≈15-cm total length (TL) mainly because of the relatively large mesh of the trawl net and possibly also because of the exclusion of areas <20-m deep from the survey area (Yeung and Yang, 2017). Despite these limitations, the bottom trawl survey offers the best available data and the only continuous time series for investigating juvenile flatfish. For this study, fish of ≤15-cm TL were classified as juveniles. According to the survey’s length-age database, NRS of ≤15-cm TL were 94% age ≤4 (16% age 2, 51% age 3, and 26% age 4), and YFS were 92% age ≤6 (20% age 3, 35% age 4, 26% age 5, and 10% age 6). Within this length class, 91% of NRS and 93% of YFS were of ≥10-cm TL. For brevity, hereinafter the terms NRS and YFS implicitly refer to juveniles of the species, unless specified otherwise (e.g. adult, all ages). Biological indices The catch per unit effort (density, no. ha−1) by station from 1982 to 2017 was used to map the spatial distributions of NRS and YFS and calculate the biological indices of spatial distribution: area of 95% presence (area), longitude (cgx) and latitude (cgy) of the centre of distribution, global (gic), and local (lic) colocation indices (Table 1;Calenge, 2006; Renard et al., 2016; Petitgas et al., 2017). The annual weight of total catch of NRS and YFS of all ages and the annual estimated number of juveniles (age ≤4) and numbers at each age from 5 to 10 were used as biological indices of abundance (Table 1). The percentage of age-1 in the juvenile class was negligible for NRS and practically zero for YFS. The mean percentage composition of ages 2, 3, and 4 in the total number of juveniles were 7, 36, and 57% for NRS and 1, 17, and 82% for YFS. Environmental indices Bottom temperature, sea ice, and wind stress indices were used to characterize annual EBS environment (Table 1). The ice cover anomaly index (ICI), ice retreat index (IRI), spring and winter wind stress anomaly indices (WS_MJ, WS_NA) were obtained from the Bering Sea climate data website (https://www.beringclimate.noaa.gov) maintained by the NOAA Pacific Marine Environmental Laboratory. The mean temperature of EBS northern inner shelf stratum 20 (BT20) and southern middle shelf stratum 31 (BT31) were obtained from the bottom trawl survey (Figure 1). Exploratory analysis had shown that BT20 and BT31 were highly correlated with bottom temperature in other EBS strata, with the EBS shelf as a whole, and with the area of the cold pool. Data analysis The influence of the EBS environment on the spatial distribution and abundance of NRS and YFS and the link between juvenile and adult abundances were examined using time-series correlation and multiple linear regressions. All analyses were performed in the R environment (R Core Team, 2018). Time series correlations The relationships among and between environmental and biological indices were examined with correlation analysis (Pearson r). Biological indices (spatial distribution and abundance, Table 1) were cross-correlated with environmental indices at a maximum lag of 4 years, corresponding with the maximum age of the juvenile class. The total number (abundance) of juveniles estimated from the survey was cross-correlated with the number at each age from 5 to 10, and with the total survey catch of all ages. The maximum lag examined for this set of correlations was 6 years, corresponding with the difference between the maximum age of juveniles and the maximum age of adults considered in this study. The modified Chelton method was used to adjust for autocorrelation in time series (Pyper and Peterman, 1998). Significant correlations (p < 0.05) were further rated for strength (5—decisive; 4—very strong; 3—strong; 2—substantial; 1—anecdotal) using the Bayes Factor (BF) in a Bayesian hypothesis test (Wetzels and Wagenmakers, 2012; Nuijten et al., 2015), with BF = 3–5 collectively classified in this study as “strong.” Species–environment models All-possible-subsets (without interaction) ordinary least squares (OLS) regression (Hebbali, 2018) was applied for each biological index (response variable) using the environmental indices (covariates) with which it was significantly (p < 0.05) and strongly (BF ≥ 3) correlated. Observations up to year 2017 were used to build the model, and the model was used to predict the biological response for 2018 and possibly beyond, depending on the availability of the environmental indices for the specific year(s) required in the model. Models within all the possible subsets that had Bayesian Information Criterion (BIC) (Schwarz, 1978) less than the minimum BIC + 2 (i.e. ΔBIC < 2) are all highly supported by the data (Burnham et al., 2002). To arrive at the “best” model, we first (1) selected the model among the ΔBIC < 2 subset with the fewest covariates, then we (2) eliminated the covariates from (1) that had variance inflation factor (vif) > 4 (high multicollinearity) (Hair et al., 2010) and coefficient estimates of p(|t|) > 0.05 (non-significant effect). Lastly, the Durbin-Watson (D-W) test for first-order autocorrelation was used to test model (2) for violation of the assumptions of normal and independent errors. If autocorrelation in the residuals was significant (D-W statistic, p < 0.05), model (2) was re-estimated by generalized least squares (GLS) with first-order autoregressive AR(1) error (Fox and Weisburg, 2010), and if the ΔBIC > 2, the GLS model was adopted as the “best” model, otherwise the OLS model (2) remained the “best” model. The coefficient of determination (r2), adjusted for the number of covariates in the model, was used as a standard measure of OLS goodness-of-fit (Neter et al., 1989). An r2 value is not a standard result of GLS, but here it was estimated by piecewise structural equation modelling (Lefcheck, 2016) for an approximate comparison with OLS. Results Variability in spatial distribution and abundance The proportion of stations surveyed where NRS were present was higher than average for an extended period in 1986–1993 and in 2004–2009, peaking in 1991 (Supplementary Figure S1). In 2010–2015, the presence of NRS in the survey area was markedly below average. For YFS, there were no extended periods of anomaly, but oscillations about the average of 1- to 3-year periods. Comparing the overall presence of juveniles in the combined EBS and NBS between 2010 and 2017 (Supplementary Figure S2), NRS were present in the NBS in 2017 but not in 2010. Their overall proportional presence in 2017 was not much greater than average (Supplementary Figure S1), indicating a change in the shape of the distribution rather than an increase in the area of distribution (range). YFS was found in the NBS in both years (Supplementary Figure S2), but there was a peak in presence in 2010 that was not present in 2017 (Supplementary Figure S1), indicating that the range across the combined EBS and NBS was much wider in 2010 than in 2017. Considering the EBS only, the centre of distribution of NRS began to shift alongshelf towards the southeast in 2009 and reached furthest south in 2010–2014 (Figure 2). As the centre moved southward the range contracted (Figure 3). After 2014, the shift reversed back towards the northwest and the range expanded, including a noticeable presence near the Pribilof Islands in 2014–2015. By 2017, the centre was again at among some of the northernmost locations ever recorded, and the distribution pattern resembled that in 2008. In contrast, the spatial distribution of YFS has been relatively constant over time (Figure 2), with the centre of distribution located around 59°N and little change in the range (Figure 3). The centre was consistently further inshore than NRS and north of 58°N. The centre of distribution of YFS was furthest south in 1985 and 2002. The centres of NRS distribution were practically the same in 2010 whether NBS was included or not, whereas in 2017 the centre for the combined EBS and NBS was north and inshore of the centre for the EBS only (Figure 2). The centres of YFS distribution for the combined EBS and NBS in 2010 and 2017 were distinctly further north and inshore of all the centres for EBS only over the years. Figure 2. Open in new tabDownload slide The centres of distribution of NRS and YFS juveniles from 1982 to 2017 for the EBS survey area only, marked by the last two digits of the year (NRS in italics; YFS in regular font). For YFS, a polygon encloses the centres of all years except those that can be distinctly separated and labelled. In addition, centres with the NBS survey area included in 2010 and 2017 are also shown as underlined two-digit year (10 and 17). For YFS, the 2010 and 2017 centres for the combined EBS and NBS are distinctly north and inshore of the 2010 and 2017 centres for EBS only, which are within the polygon. For NRS, the centres for the combined EBS and NBS and for EBS only coincide in 2010, whereas the centre for the combined EBS and NBS is north and inshore of the centre for EBS only in 2017. Figure 2. Open in new tabDownload slide The centres of distribution of NRS and YFS juveniles from 1982 to 2017 for the EBS survey area only, marked by the last two digits of the year (NRS in italics; YFS in regular font). For YFS, a polygon encloses the centres of all years except those that can be distinctly separated and labelled. In addition, centres with the NBS survey area included in 2010 and 2017 are also shown as underlined two-digit year (10 and 17). For YFS, the 2010 and 2017 centres for the combined EBS and NBS are distinctly north and inshore of the 2010 and 2017 centres for EBS only, which are within the polygon. For NRS, the centres for the combined EBS and NBS and for EBS only coincide in 2010, whereas the centre for the combined EBS and NBS is north and inshore of the centre for EBS only in 2017. Figure 3. Open in new tabDownload slide Interannual variability in the spatial distributions of NRS and YFS juveniles in the EBS from 2008 to 2017. The centre of the ellipse delineates the mean location, and the major and minor axes represent the maximum and minimum inertia, respectively. The area of the ellipse is proportional to the total inertia. Symbols are centred at survey stations and sized proportionally to the relative density for the year. Figure 3. Open in new tabDownload slide Interannual variability in the spatial distributions of NRS and YFS juveniles in the EBS from 2008 to 2017. The centre of the ellipse delineates the mean location, and the major and minor axes represent the maximum and minimum inertia, respectively. The area of the ellipse is proportional to the total inertia. Symbols are centred at survey stations and sized proportionally to the relative density for the year. The range of NRS on the EBS shelf was twice that of YFS. The range of both NRS and YFS generally fluctuated about the average with no distinct multi-year patterns, except for NRS there was one extended depression in their range in 2008–2013 (Supplementary Figure S3). The spatial overlap of NRS and YFS distributions were high over the EBS inner shelf in general (Figure 3) as reflected in the high (0.8 ± 0.2; mean ± 1 SD) global index of colocation (gic), but the low (0.3 ± 0.2) local index of colocation (lic) suggested that within the inner shelf NRS and YFS co-occurred at the same locations only periodically, in an ∼3- to 6-year cycle, with a peak in 1985 (Supplementary Figure S4). The variability in colocation was driven by the variability in the spatial distribution of NRS. The two populations were most distinctly separated for an extended period between 2009 and 2015 (Figure 3). Correlations among biological and environmental indices The expansion/contraction of the distribution range and the shift in the centre of distribution were tightly coupled for both NRS and YFS. The contraction of the range of NRS led the shift in the centre of distribution to the southeast by 0–3 years (r = 0.56–0.79, p = 0.0001–0.016). The expansion of their range was correlated with warming bottom temperature over the EBS shelf (r = 0.50, p = 0.01). The range expansion of YFS led the eastward shift in the centre of distribution by 2 years (r = 0.52, p < 0.001), but was not correlated with bottom temperature over the EBS shelf. The mean summer bottom temperature of the entire EBS shelf was significantly correlated with its component inner and middle shelf strata (r = 0.75–0.91, p < 0.002). Below average EBS shelf bottom temperature (2006–2010) (Figure 4a) led the shift in the centre of distribution of NRS eastward (inshore) and northward by 1–3 years (r = 0.45–0.56, p = 0.005–0.031) (Figure 4b and c). The correlation was highest when temperature led by 2 years. Higher temperature also significantly led higher number of NRS by 1–3 years (r = 0.43–0.48, p < 0.005). The estimated number of NRS went through larger magnitudes and longer periods of anomalies about the 1982–2017 average since 1997 (Figure 4d). The latest negative anomaly between 2008 and 2015 corresponded with the southeast shift in the population. The number of YFS was much lower than that of NRS, and anomalies were not as large. Their abundance was below average for most of the time series, with occasional spikes of positive anomaly, notably in 1982, 2007, and 2017 (Figure 4d). Mean temperature across the EBS shelf was not significantly correlated with the centre of distribution nor the number of YFS (Figure 4a–d). Figure 4. Open in new tabDownload slide (a) The mean bottom temperature on the EBS shelf, (b) the longitude and (c) latitude of the centres of distribution and (d) the estimated number of NRS and YFS juveniles from 1982 to 2017. The horizontal line depicts the series mean. Also shown for 2010 and 2017 in (b–d) are the values with the NBS area included (hollow markers). Figure 4. Open in new tabDownload slide (a) The mean bottom temperature on the EBS shelf, (b) the longitude and (c) latitude of the centres of distribution and (d) the estimated number of NRS and YFS juveniles from 1982 to 2017. The horizontal line depicts the series mean. Also shown for 2010 and 2017 in (b–d) are the values with the NBS area included (hollow markers). None of the biological indices for YFS was significantly cross-correlated with any of the environmental indices at a lag of up to 4 years. In contrast, most of the biological indices for NRS were strongly correlated with environmental indices at a lag of up to 3 years, and correlations with the ice indices were particularly strong (Supplementary Table S1a). However, wind stress indices were not significantly correlated with any of the biological indices for NRS. The centre of distribution (longitude and latitude) of NRS was strongly correlated with the ICI, the IRI, and the bottom temperature in stratum 31 (BT31) of the prior 3 years. High and prolonged ice cover and low bottom temperature led the shift of the population centre eastward (inshore) and southward. High bottom temperature led the shift of the population centre northward. The area of 95% presence (area) decreased with increased ice and decreased temperature. The gic had an overall stronger correlation with environmental indices of the prior 3 years than the lic. Overlap of the distributions of NRS and YFS decreased with increased ice cover and decreased temperature. The number of NRS (NRS1-4) was strongly correlated with the ICI of 2–3 years ago and the BT31 of 3 years ago (Supplementary Table S1a). Low ice cover and high bottom temperature led large number of NRS. The number at each adult age from 5 to 9 (NRS5–NRS9) was correlated with NRS1-4 at the expected lags within an error in age estimation of ±1 year. For example, NRS5 was strongly correlated with NRS1-4 at a 0- to 2-year lag (Supplementary Table S1b). The correlation between NRS1-4 and NRS10 was not significant. The number of YFS2-4 was not correlated with the number of YFS at any age ≥5. The catch of NRS significantly lagged age-5 (by 3 years) to age-8 (by 0 years) numbers. The magnitude of the correlation coefficient r for the specific lag at age suggested that the catch was mainly composed of age-8 (Supplementary Table S1c). The catch of YFS significantly lagged age-5 (by 2 years) to age-10 (by 0 years). Similar to NRS, the catch of YFS appeared to be mainly composed of ages 7–8, but included a wider range of ages. Ice and bottom temperature indices were significantly and strongly correlated at no time lag (rICI, IRI = 0.75, rICI, BT20 = −0.47, rICI, BT31 = −0.82, rIRI, BT20 = −0.52, rIRI, BT31 = −0.84, rBT20, BT31 = 0.58). When lag-transformed, these environmental indices each showed different degrees of correlation with biological indices (Supplementary Table S1a). Ice cover (ICI) had the strongest correlation with all the biological variables (|r| = 0.53–0.87, leading by 0–3 years). Bottom temperature on the southern middle shelf (BT31) was also strongly correlated with the distribution and abundance of NRS (|r| = 0.47–0.60, leading by 2–3 years). Bottom temperature on the northern inner shelf (BT20) was only strongly correlated with the latitude of the centre of distribution of NRS (|r| = 0.56–0.60, leading by 2–3 years). Species–environment models Bottom temperature and ice indices were strong environmental covariates for building informative regression models of biological indices for NRS (Supplementary Table S2, Figure 5). Only the global colocation index model was re-estimated with GLS, whereas all other biological indices for NRS were modelled with OLS. Figure 5. Open in new tabDownload slide Linear regression models fitted (bold line) to observed (filled diamond) juvenile NRS (a) longitude—cgx, and (b) latitude—cgy, of the centre of distribution, (c) abundance—NRS1-4, (d) area of 95% presence—area, and the indices of colocation between the juveniles of NRS and YFS: (e) local—lic and (f) global—gic, in year T using environmental covariates (ice cover index—ICI; ice retreat index—IRI; mean bottom temperature at stratum j—BTj). All models were fitted by OLS except gic, which was fitted by GLS. The 95% confidence bands (light line) are shown. The standard adjusted r2 measures the goodness-of-fit of OLS models, and an r2 for the GLS model was separately estimated by piecewise structural equation modelling. Values beyond 2017 (open circle) were predicted by the model. See Methods section and Supplementary Table S2. Figure 5. Open in new tabDownload slide Linear regression models fitted (bold line) to observed (filled diamond) juvenile NRS (a) longitude—cgx, and (b) latitude—cgy, of the centre of distribution, (c) abundance—NRS1-4, (d) area of 95% presence—area, and the indices of colocation between the juveniles of NRS and YFS: (e) local—lic and (f) global—gic, in year T using environmental covariates (ice cover index—ICI; ice retreat index—IRI; mean bottom temperature at stratum j—BTj). All models were fitted by OLS except gic, which was fitted by GLS. The 95% confidence bands (light line) are shown. The standard adjusted r2 measures the goodness-of-fit of OLS models, and an r2 for the GLS model was separately estimated by piecewise structural equation modelling. Values beyond 2017 (open circle) were predicted by the model. See Methods section and Supplementary Table S2. The ice indices were covariates in all the best models of NRS abundance and distribution. The models fit spatial distribution best (adjusted r2 > 0.7, Figure 5a and b), and abundance (Figure 5c), presence (Figure 5d) and the colocation indices (Figure 5e and f) moderately well (adjusted r2 = 0.4–0.6). Bottom temperature was significant only in the models for abundance and longitude. There was no strong indication that the biological responses of NRS predicted by models for 2018 would be much different from those observed in 2017. In 2018, the predicted NRS abundance remained high (Figure 5c), likewise the NRS distribution range and the spatial overlap between NRS and YFS—although they trended slightly downwards (Figure 5d–f), and the centre of distribution of NRS remained far to the northwest (Figure 5a and b). Discussion Ice cover and bottom temperature in the preceding 1–3 years are the best predictors of the abundance and distribution of NRS in the EBS in this study. Low ice cover and high bottom temperature characteristic of a warm stanzas are strongly correlated with higher abundance of NRS, supporting previous observations of positive correlation between temperature and NRS densities (Cooper et al., 2014; Cooper and Nichol, 2016; Wilderbuer and Nichol, 2016). High ice cover and low bottom temperature characteristic of cold stanzas are correlated with lower abundance of NRS, and with the southeastward displacement and contraction of their distribution range. The high abundance and northerly distribution of predominantly ages 3–4 NRS recently observed in the EBS would correspond with the current warm thermal stanza that began in 2014. In contrast, variability in the abundance and distribution of YFS is relatively low, and not significantly correlated with ice or temperature indices. However, opposite patterns have been observed in the adult populations (Nichol et al., 2019), where temperature indices were positively correlated with the biomass estimates of YFS adults, but had no correlation with the biomass of NRS adults, suggesting that ecological dynamics and spatial availability are specific to species and life stages. It was hypothesized that colder bottom temperature delay migration and spawning of YFS, causing mature individuals to reside in nearshore nursery grounds that are unavailable to the survey, thus decreasing the biomass estimates during cold years (Nichol et al., 2019). The low correlations between temperature and both adult NRS biomass and juvenile YFS abundance may be partly due to neither one undergoing significant cross-shore migration in and out of the survey area during the survey period. Generally, juveniles of NRS and YFS are known to be distributed on the EBS inner shelf (Bartolino et al., 2011), but this study finds that the distribution of YFS is centred inshore of NRS. Higher spatial overlap of the two species during warm stanzas resulted from the spread of NRS northward along the inner shelf; lower overlap resulted from the contraction of NRS distribution southward during cold stanzas. Many temperature-related factors during the early life history of flatfish, such as the timing of spawning, larval duration and transport, nursery conditions, etc., can affect spatial distribution of juvenile flatfish (Teal et al., 2012; Duffy-Anderson et al., 2015; Laurel et al., 2015; van der Veer et al., 2015). Juvenile YFS may be less affected by the fluctuations in the EBS thermal environment than NRS because of their different early life histories. NRS settle into nearshore habitat in late spring or early summer (Norcross et al., 1995), when the bottom temperature tends to be variable (Overland et al., 1999), and during cold stanzas the influence of the cold pool over this area and during this time of the year is still relatively strong (Stabeno et al., 2012). When the cold pool is prominent and expansive, it could limit the offshore and northern extent of suitable thermal environment for settlement. YFS settle later in summer (Wilderbuer et al., 1992), when the influence of the cold pool is weaker and bottom temperature nearshore is relatively stable and high (Bakkala, 1981; Yeung and Yang, 2018). NRS may also have a greater dependence on larval transport for successful settlement. Putative spawning locations of NRS occur along the Alaska Peninsula (Lanksbury et al., 2007) distant from their nursery habitat, so their larvae are dependent on current transport, which varies between warm and cold stanzas (Stabeno et al., 2012). In contrast, YFS migrate across the shelf to spawn near their nursery habitat (Nichol and Acuna, 2001; Norcross and Holladay, 2005), so their larvae may be less susceptible to variable transport. Springtime winds and drift patterns, because of their correlation with surface currents, are used as indicators of recruitment for winter-spawning flatfish such as NRS (Wilderbuer and Nichol, 2016; North Pacific Fishery Management Council, 2017b). In this study, however, wind stress indices are not significantly correlated with NRS abundance and distribution. According to the models in this study, the present and continuing warm stanza is expected to produce high NRS abundance up to 3 years past its dissolution, and the range is expected to remain extended towards the north up to 1 year after. The NRS fishery is composed of fish aged primarily 6–10. As NRS abundance is strongly and positively correlated with the abundance of subsequent age classes up to age 9, the adult biomass is also expected to be high in the next few years as the age 3–4 fish from 2016 to 2018 recruit into the fishery. There is no significant relationship between YFS abundance and distribution and the thermal environment. The true relationship could have been obscured by inadequate sampling coverage of their nearshore habitat (Nichol, 1997; Nichol et al., 2019). It is less likely to have been smoothed over by the large number of year classes in the defined juvenile size range of ≤15-cm TL (YFS age ≤6; NRS age ≤4), as the models in this study showed no significant relationships even when the juvenile class basically consisted of only ages 3–4. Although juvenile YFS abundance has been increasing since 2015, the future population trajectory is unpredictable, as juvenile YFS abundance is not correlated with the abundance of subsequent adult age classes, nor does it correspond well with stock assessment biomass estimates for the total population (Wilderbuer et al., 2017; Nichol et al., 2019). Whereas juvenile YFS abundance has been below average in 1983–2003, and 2013–2015, total YFS biomass shows a mostly upward trend since 1990, except for a downturn in 2012—when the second coldest summer bottom temperature ever was recorded on the shelf, and in 2015. The biomass peaked in 2016 during the highest bottom temperature ever recorded. Since 2000, biomass estimates have been positively correlated with shelf bottom temperature. This correlation is perhaps more a reflection of the availability of adult YFS to the survey area than an increase in population abundance, as the migration of adult YFS from the shelf break to shallower waters in summer better matched the timing and area of the survey (Wilderbuer et al., 2017; Nichol et al., 2019). Warmer Bering Sea temperatures and a diminished cold pool led the northward expansion of the distribution range of NRS juveniles. The northward trend was evident in 2014 at the start of the present warm stanza, and confirmed by the presence of NRS in the NBS in 2017. They were not caught in the NBS in 2010 (Stevenson and Lauth, 2019). The 2017 juvenile class was the product of a warm stanza, whereas the 2010 juvenile class was the product of a cold stanza and apparently did not reach as far north. The NBS surveys in 2010 and 2017 both showed the presence of YFS, suggesting that they may be endemic to the NBS inner shelf, which experiences high summer bottom temperatures and has rich prey resources (Hamazaki et al., 2005; Yeung and Yang, 2014). Considering the entire US Bering Sea shelf, the proportions of YFS in the NBS were small in 2010 and 2017, and its contribution to the overall production of the YFS population is unknown. To evaluate how changing thermal environment would affect the abundance and distribution of NRS and YFS, not only each species’ own physiological responses to temperature, but the responses of other linked species in the ecosystem need to be considered. The competitors and predators of flatfishes may also spread northward (van Hal et al., 2010), and with increased growth and consumption rates (Akimova et al., 2016). Warmer environment increases the spatial overlap between the ecologically similar NRS and YFS, potentially increasing competition. The biomass of Pacific cod (Gadus macrocephalus) and pollock, two of the main predators of the juvenile flatfishes, respectively increased 900 and 6000% in the NBS in 2017 since the previous survey there in 2010 (Stevenson and Lauth, 2019), while the population of Pacific halibut (Hippoglossoides stenolepis), whose main prey are YFS, declined by more than 20% (Lauth, 2011; Lauth et al., 2019). The timing and extent of ice cover and the cold pool affect phytoplankton and zooplankton production (Overland and Stabeno, 2004), and probably also the production of the macrobenthos on which the flatfishes prey (Lin et al., 2019). Responses observed thus far in the managed EBS groundfish species were not only variable between feeding guilds, but also within each type of stanza. For example, in the present warm stanza, the productivity of apex predators such as arrowtooth flounder (Atheresthes stomias) increased, whereas pelagic foragers such as pollock (Gadus chalcogrammus) was neutral, differing from the previous warm stanza when the productivity of both species decreased. Among benthic foragers, the productivity of NRS declined 25% in 2015 after the recent cold stanza, but YFS increased 50% in 2016 (North Pacific Fishery Management Council, 2018). Until there is sufficient understanding of the underlying ecological mechanisms, the responses of groundfish species to alternating thermal stanzas and to the warming climate will not be able to be predicted with confidence (Mueter and Litzow, 2008; Hunt et al., 2011; Teal et al., 2012; Duffy-Anderson et al., 2017). Acknowledgements We thank all the collectors and custodians of the AFSC bottom trawl survey data and especially Dan Nichol for help with extracting the flatfish data for this study. References Akimova A. , Hufnagl M. , Kreus M. , Peck M. A. 2016 . Modeling the effects of temperature on the survival and growth of North Sea cod (Gadus morhua) through the first year of life . Fisheries Oceanography , 25 : 193 – 209 . Google Scholar Crossref Search ADS WorldCat Armistead C. , Nichol D. G. 1993 . 1990 bottom trawl survey of the eastern Bering Sea and continental shelf. NTIS No. PB93-156677: 159 pp. Bakkala R. G. 1981 . Population characteristics and ecology of yellowfin sole. In The Eastern Bering Sea Shelf: Oceanography and Resources , pp. 553 – 574 . Ed. by Hood D. W. , Calder J. A. . University of Washington Press , Seattle . Google Preview WorldCat COPAC Bartolino V. , Ciannelli L. , Bacheler N. M. , Chan K.-S. 2011 . Ontogenetic and sex-specific differences in density-dependent habitat selection of a marine fish population . Ecology , 92 : 189 – 200 . Google Scholar Crossref Search ADS PubMed WorldCat Benson A. J. , Trites A. W. 2002 . Ecological effects of regime shifts in the Bering Sea and eastern North Pacific Ocean . Fish and Fisheries , 3 : 95 – 113 . Google Scholar Crossref Search ADS WorldCat Bohle M. S. , Bakkala R. G. 1994 . Data Report: 1978 bottom trawl survey of the eastern Bering Sea groundfish. U.S. Department of Commerce, NOAA Technical Memorandum, NMFS F/NWC-55: 164 pp. Britten G. L. , Dowd M. , Worm B. 2016 . Changing recruitment capacity in global fish stocks . Proceedings of the National Academy of Sciences of the United States of America , 113 : 134 – 139 . Google Scholar Crossref Search ADS PubMed WorldCat Burnham K. P. , Anderson D. R. , Burnham K. P. 2002 . Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach . Springer , New York . Google Preview WorldCat COPAC Calenge C. 2006 . The package “adehabitat” for the R software: a tool for the analysis of space and habitat use by animals . Ecological Modelling , 197 : 516 – 519 . Google Scholar Crossref Search ADS WorldCat Coachman L. K. 1986 . Circulation, water masses, and fluxes on the southeastern Bering Sea shelf . Continental Shelf Research , 5 : 23 – 108 . Google Scholar Crossref Search ADS WorldCat Conner J. , Lauth R. R. 2017 . Results of the 2016 eastern Bering Sea continental shelf bottom trawl survey of groundfish and invertebrate resources. U.S. Department of Commerce. NOAA Technical Memorandum, NMFS-AFSC-352. 159, pp. Cooper D. W. , Duffy-Anderson J. T. , Norcross B. L. , Holladay B. A. , Stabeno P. 2014 . Northern rock sole (Lepidopsetta polyxystra) juvenile nursery areas in the eastern Bering Sea in relation to hydrography and thermal regimes . ICES Journal of Marine Science , 71 : 1683 – 1695 . Google Scholar Crossref Search ADS WorldCat Cooper D. W. , Duffy-Anderson J. T. , Stockhausen W. T. , Cheng W. 2013 . Modeled connectivity between northern rock sole (Lepidopsetta polyxystra) spawning and nursery areas in the eastern Bering Sea . Journal of Sea Research , 84 : 2 – 12 . Google Scholar Crossref Search ADS WorldCat Cooper D. W. , Nichol D. G. 2016 . Juvenile northern rock sole (Lepidopsetta polyxystra) spatial distribution and abundance patterns in the eastern Bering Sea: spatially-dependent production linked to temperature . ICES Journal of Marine Science , 73 : 1138 – 1146 . Google Scholar Crossref Search ADS WorldCat Duffy-Anderson J. T. , Bailey K. M. , Cabral H. N. , Nakata H. , Van der Veer H. W. 2015 . The planktonic stages of flatfishes: physical and biological interactions in transport processes. In Flatfishes: Biology and Exploitation , 2nd edn, pp. 132 – 170 . Ed. by Gibson R. N. , Nash R. D. M. , Geffen A. J. , Van der Veer H. W. . John Wilery & Sons, Ltd ., West Sussex, UK . Google Preview WorldCat COPAC Duffy-Anderson J. T. , Stabeno P. J. , Siddon E. C. , Andrews A. G. , Cooper D. W. , Eisner L. B. , Farley E. V. et al. 2017 . Return of warm conditions in the southeastern Bering Sea: phytoplankton - Fish . PLoS One , 12 : e0178955. Google Scholar Crossref Search ADS PubMed WorldCat Fedewa E. J. , Miller J. A. , Hurst T. P. , Jiang D. 2017 . The potential effects of pre-settlement processes on post-settlement growth and survival of juvenile northern rock sole (Lepidopsetta polyxystra) in Gulf of Alaska nursery habitats . Estuarine Coastal and Shelf Science , 189 : 46 – 57 . Google Scholar Crossref Search ADS WorldCat Fossheim M. , Primicerio R. , Johannesen E. , Ingvaldsen R. B. , Aschan M. M. , Dolgov A. V. 2015 . Recent warming leads to a rapid borealization of fish communities in the Arctic . Nature Climate Change , 5 : 673. Google Scholar Crossref Search ADS WorldCat Fox J. , Weisberg S. 2010 . Time-series Regression and Generalized Least Squares: An Appendix to An R Companion to Applied Regression, 2nd edn. https://socialsciences.mcmaster.ca/jfox/Books/Companion-2E/appendix/Appendix-Timeseries-Regression.pdf (last accessed 20 September 2019). Grebmeier J. M. , Overland J. E. , Moore S. E. , Farley E. V. , Carmack E. C. , Cooper L. W. , Frey K. E. et al. 2006 . A major ecosystem shift in the northern Bering Sea . Science , 311 : 1461 – 1464 . Google Scholar Crossref Search ADS PubMed WorldCat Hair J. , Black W. C. , Babin B. J. , Anderson R. E. 2010 . Multivariate data analysis. Pearson Education International , Upper Saddle River , New Jersey . Google Preview WorldCat COPAC Hamazaki T. , Fair L. , Watson L. , Brennan E. 2005 . Analyses of Bering Sea bottom-trawl surveys in Norton Sound: absence of regime shift effect on epifauna and demersal fish . ICES Journal of Marine Science , 62 : 1597 – 1602 . Google Scholar Crossref Search ADS WorldCat Hebbali A. 2018 . olsrr: Tools for building OLS regression models. R package version 0.5.1. https://CRAN.R-project.org/package=olsrr Hollowed A. B. , Angliss R. P. , Sigler M. F. , Megrey B. A. , Ito D. H. 2007 . Implementation plan for Loss of Sea Ice (LOSI) program. AFSC Processed Rep. 2007–5, 48 p. Alaska Fish. Sci. Cent., NOAA, Natl. Mar, Fish. Serv., Seattle, WA. Hollowed A. B. , Planque B. , Loeng H. 2013 . Potential movement of fish and shellfish stocks from the sub-Arctic to the Arctic Ocean . Fisheries Oceanography , 22 : 355 – 370 . Google Scholar Crossref Search ADS WorldCat Hunt G. L. , Coyle K. O. , Eisner L. B. , Farley E. V. , Heintz R. A. , Mueter F. , Napp J. M. et al. 2011 . Climate impacts on eastern Bering Sea foodwebs: a synthesis of new data and an assessment of the Oscillating Control Hypothesis . ICES Journal of Marine Science , 68 : 1230 – 1243 . Google Scholar Crossref Search ADS WorldCat Intergovernmental Panel on Climate Change. 2014 . Anthropogenic and natural radiative forcing. In Climate Change 2013 – the Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change , pp. 659 – 740 . Cambridge University Press , Cambridge . WorldCat COPAC Kinder T. H. , Schumacher J. D. 1981 . Hydrographic structure over the continental shelf of the southeastern Bering Sea. In The Eastern Bering Sea Shelf: Oceanography and Resources , pp. 31 – 52 . Ed. by Hood D. W. , Calder J. A. . University of Washington Press , Seattle . Google Preview WorldCat COPAC Lanksbury J. A. , Duffy-Anderson J. T. , Mier K. L. , Busby M. S. , Stabeno P. J. 2007 . Distribution and transport patterns of northern rock sole, Lepidopsetta polyxystra, larvae in the southeastern Bering Sea . Progress in Oceanography , 72 : 39 – 62 . Google Scholar Crossref Search ADS WorldCat Laurel B. , Basilio A. , Danley C. , Ryer C. , Spencer M. 2015 . Substrate preference and delayed settlement in northern rock sole larvae Lepidopsetta polyxystra . Marine Ecology Progress Series , 519 : 183 – 193 . Google Scholar Crossref Search ADS WorldCat Laurel B. J. , Danley C. , Haines S. 2014 . The effects of temperature on growth, development and settlement of northern rock sole larvae (Lepidopsetta polyxystra) . Fisheries Oceanography , 23 : 495 – 505 . Google Scholar Crossref Search ADS WorldCat Lauth R. R. 2011 . Results of the 2010 eastern and northern Bering Sea continental shelf bottom trawl survey of groundfish and invertebrate fauna. U.S. Department of Commerce. NOAA Technical Memorandum, NMFS-AFSC-227: 256, pp. Lauth R. R. , Dawson E. J. , Conner J. 2019 . Results of the 2017 eastern and northern Bering Sea continental shelf bottom trawl survey of groundfish and invertebrate fauna. U.S. Department of Commerce. NOAA Technical Memorandum, NMFS-AFSC-396: 260, pp. Lefcheck J. S. 2016 . piecewiseSEM: piecewise structural equation modelling in r for ecology, evolution, and systematics . Methods in Ecology and Evolution , 7 : 573 – 579 . Google Scholar Crossref Search ADS WorldCat Lin H. , Liu K. , Wang J. , Lin J. , Huang Y. , He X. , Li Z. et al. 2019 . Spatial variability of macrobenthic production in the Bering Sea . Polar Biology , 42 : 449 – 460 . Google Scholar Crossref Search ADS WorldCat Link J. S. 2016 . The importance of understanding ecosystem processes in the Eastern Bering Sea . Deep-Sea Research Part II: Topical Studies in Oceanography , 134 : 424 – 425 . Google Scholar Crossref Search ADS WorldCat Matarese A. C. , Blood D. M. , Picquelle S. J. , Benson J. L. 2003 . Atlas of abundance and distribution patterns of ichthyoplankton from the northeast Pacific Ocean and Bering Sea ecosystems based on researched conducted by the Alaska Fisheries Science Center (1972-1996). NOAA Prof. Paper NMFS, 1: 281. pp. Mueter F. J. , Litzow M. A. 2008 . Sea ice retreat alters the biogeography of the Bering Sea continental shelf . Ecological Applications , 18 : 309 – 320 . Google Scholar Crossref Search ADS PubMed WorldCat Neter J. , Wasserman W. , Kutner M. H. 1989 . Applied Linear Regression Models . Irwin Press , Homewood, IL . 688 pp. Google Preview WorldCat COPAC Nichol D. G. 1997 . Effects of geography and bathymetry on growth and maturity of yellowfin sole, Pleuronectes asper, in the eastern Bering Sea . Fishery Bulletin , 95 : 494 – 503 . WorldCat Nichol D. G. , Acuna E. I. 2001 . Annual and batch fecundities of yellowfin sole, Limanda aspera, in the eastern Bering Sea . Fishery Bulletin , 99 : 108 – 122 . WorldCat Nichol D. G. , Kotwicki S. , Wilderbuer T. K. , Lauth R. R. , Ianelli J. N. 2019 . Availability of yellowfin sole Limanda aspera to the eastern Bering Sea trawl survey and its effect on estimates of survey biomass . Fisheries Research , 211 : 319 – 330 . Google Scholar Crossref Search ADS WorldCat Norcross B. L. , Holladay B. A. 2005 . Feasibility to design and implement a nearshore juvenile flatfish survey – eastern Bering Sea. Final Technical Report to the Cooperative Institute for Arctic Research Award #NA17RJ1224. Norcross B. L. , Holladay B. A. , Müter F. J. 1995 . Nursery area characteristics of pleuronectids in coastal Alaska, USA . Netherlands Journal of Sea Research , 34 : 161 – 175 . Google Scholar Crossref Search ADS WorldCat North Pacific Fishery Management Council. 2017a . 2017 North Pacific groundfish stock assessment and fishery evaluation reports for 2018 fisheries. https://www.fisheries.noaa.gov/alaska/population-assessments/2017-north-pacific-groundfish-stock-assessments (last accessed 20 September 2019). North Pacific Fishery Management Council. 2017b . Ecosystem Considerations 2017. Status of the Eastern Bering Sea Marine Ecosystem. https://repository.library.noaa.gov/view/noaa/19464 (last accessed 20 September 2019). North Pacific Fishery Management Council. 2018. Ecosystem Status Report 2018. Eastern Bering Sea. December 2018. https://access.afsc.noaa.gov/reem/ecoweb/pdf/2018ecosysEBS-508.pdf (last accessed 21 September 2019). Nuijten M. B. , Wetzels R. , Matzke D. , Dolan C. V. , Wagenmakers E.-J. 2015 . BayesMed: default Bayesian hypothesis tests for correlation, partial correlation, and mediation . R package version 1.0.1. https://CRAN.R-project.org/package=BayesMed Overland J. E. , Alheit J. , Bakun A. , Hurrell J. W. , Mackas D. L. , Miller A. J. 2010 . Climate controls on marine ecosystems and fish populations . Journal of Marine Systems , 79 : 305 – 315 . Google Scholar Crossref Search ADS WorldCat Overland J. E. , Salo S. A. , Kantha L. H. , Clayson C. A. 1999 . Thermal stratification and mixing on the Bering Sea shelf. In Dynamics of the Bering Sea , pp. 129 – 216 . Ed. by Loughlin T. R. , Ohtani K. . University of Alaska Sea Grant, AK-SG-99-03 , Fairbanks, Alaska , 838 pp. Google Preview WorldCat COPAC Overland J. E. , Stabeno P. J. 2004 . Is the climate of the Bering Sea warming and affecting the ecosystem? Eos , 85 : 309 – 316 . Google Scholar Crossref Search ADS WorldCat Petitgas P. , Woillez M. , Rivoirard J. , Renard D. , Bez N. 2017 . Handbook of geostatistics in R for fisheries and marine ecology. ICES Cooperative Research Report, No. 177. Pyper B. J. , Peterman R. M. 1998 . Comparison of methods to account for autocorrelation in correlation analyses of fish data . Canadian Journal of Fisheries and Aquatic Sciences , 55 : 2127 – 2140 . Google Scholar Crossref Search ADS WorldCat R Core Team. 2018 . R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing , Vienna, Austria . WorldCat COPAC Rätz H.-J. , Lloret J. 2003 . Variation in fish condition between Atlantic cod (Gadus morhua) stocks, the effect on their productivity and management implications . Fisheries Research , 60 : 369 – 380 . Google Scholar Crossref Search ADS WorldCat Renard D. , Bez N. , Desassis N. , Beucher H. , Ors F. , Freulon X. 2016 . RGeostats: Geostatistical package. R package version 11.0.5. http://cg.ensmp.fr/rgeostats Schwarz G. 1978 . Estimating the dimension of a model . Annals of Statistics , 6 : 461 – 464 . Google Scholar Crossref Search ADS WorldCat Sigler M. F. , Renner M. , Danielson S. L. , Eisner L. B. , Lauth R. R. , Kuletz K. J. , Logerwell E. A. et al. 2011 . Fluxes, fins, and feathers: relationships among the Bering, Chukchi, and Beaufort Seas in a time of climate change . Oceanography , 24 : 250 – 265 . Google Scholar Crossref Search ADS WorldCat Stabeno P. , Duffy-Anderson J. , Eisner L. , Farley E. , Heintz R. , Mordy C. 2017 . Return of warm conditions in the southeastern Bering Sea: physics to fluorescence . PLoS One , 12 : e0185464 . Google Scholar Crossref Search ADS PubMed WorldCat Stabeno P. J. , Kachel N. B. , Moore S. E. , Napp J. M. , Sigler M. , Yamaguchi A. , Zerbini A. N. 2012 . Comparison of warm and cold years on the southeastern Bering Sea shelf and some implications for the ecosystem. Deep-Sea Research Part II: Topical Studies in Oceanography, 65–70: 31–45. Stark J. W. 2012 . Contrasting maturation and growth of northern rock sole in the eastern Bering Sea and Gulf of Alaska for the purpose of stock management . North American Journal of Fisheries Management , 32 : 93 – 99 . Google Scholar Crossref Search ADS WorldCat Stevenson D. E. , Lauth R. R. 2019 . Bottom trawl surveys in the northern Bering Sea indicate recent shifts in the distribution of marine species . Polar Biology , 42 : 407 – 421 . Google Scholar Crossref Search ADS WorldCat Teal L. R. , van Hal R. , van Kooten T. , Ruardij P. , Rijnsdorp A. D. 2012 . Bio-energetics underpins the spatial response of North Sea plaice (Pleuronectes platessa L.) and sole (Solea solea L.) to climate change . Global Change Biology , 18 : 3291 – 3305 . Google Scholar Crossref Search ADS WorldCat van der Veer H. W. , Freitas V. , Leggett W. C. 2015 . Recruitment level and variability. In Flatfishes: Biology and Exploitation , 2nd edn, pp. 185 – 206 . Ed. by Gibson R. N. , Nash R. D. M. , Geffen A. J. , Van der Veer H. W. . John Wilery & Sons, Ltd ., West Sussex, UK . Google Preview WorldCat COPAC van der Veer H. W. , Witte J. I. J. 1993 . The ‘maximum growth/optimal food condition’ hypothesis: a test for 0-group plaice Pleuronectes platessa in the Dutch Wadden Sea . Marine Ecology Progress Series , 101 : 81 – 90 . Google Scholar Crossref Search ADS WorldCat van Hal R. , Smits K. , Rijnsdorp A. D. 2010 . How climate warming impacts the distribution and abundance of two small flatfish species in the North Sea . Journal of Sea Research , 64 : 76 – 84 . Google Scholar Crossref Search ADS WorldCat Wakabayashi K. 1989 . Studies on the fishery biology of yellowfin sole in the eastern Bering Sea . Bulletin of the Far Seas Fisheries Research Laboratory , 26 : 21 – 152 . WorldCat Wetzels R. , Wagenmakers E.-J. 2012 . A default Bayesian hypothesis test for correlations and partial correlations . Psychonomic Bulletin and Review , 19 : 1057 – 1064 . Google Scholar Crossref Search ADS PubMed WorldCat Wilderbuer T. K. , Nichol D. G. 2012 . Assessment of the northern rock sole stock in the Bering Sea and Aleutian Islands. 2012 North Pacific groundfish stock assessment and fishery evaluation reports for 2013. https://www.afsc.noaa.gov//refm/stocks/2012_assessments.htm. Wilderbuer T. K. , Nichol D. G. 2016 . Assessment of the northern rock sole stock in the Bering Sea and Aleutian Islands. 2016 North Pacific groundfish stock assessment and fishery evaluation reports for 2017. https://www.fisheries.noaa.gov/alaska/population-assessments/2016-north-pacific-groundfish-stock- assessments. Wilderbuer T. K. , Nichol D. G. , Ianelli J. 2017 . Assessment of the yellowfin sole stock in the Bering Sea and Aleutian Islands. 2017 North Pacific groundfish stock assessment and fishery evaluation reports for 2018. https://www.fisheries.noaa.gov/alaska/population-assessments/north-pacific-groundfish-stock-assessment-and-fishery-evaluation/2017-north-pacific-groundfish-stock-assessments. Wilderbuer T. K. , Walters G. E. , Bakkala R. G. 1992 . Yellowfin sole, Pleuronectes asper, of the eastern Bering Sea: biological characteristics, history of exploitation, and management . Marine Fisheries Review , 54 : 1 – 18 . WorldCat Wood K. R. , Bond N. A. , Danielson S. L. , Overland J. E. , Salo S. A. , Stabeno P. J. , Whitefield J. 2015 . A decade of environmental change in the Pacific Arctic region . Progress in Oceanography , 136 : 12 – 31 . Google Scholar Crossref Search ADS WorldCat Wyllie-Echeverria T. , Wooster W. S. 1998 . Year-to-year variations in Bering Sea ice cover and some consequences for fish distributions . Fisheries Oceanography , 7 : 159 – 170 . Google Scholar Crossref Search ADS WorldCat Yeung C. , Yang M.-S. 2014 . Habitat and infauna prey availability for flatfishes in the northern Bering Sea . Polar Biology , 37 : 1769 – 1784 . Google Scholar Crossref Search ADS WorldCat Yeung C. , Yang M.-S. 2017 . Habitat quality of the coastal southeastern Bering Sea for juvenile flatfishes from the relationships between diet, body condition and prey availability . Journal of Sea Research , 119 : 17 – 27 . Google Scholar Crossref Search ADS WorldCat Yeung C. , Yang M.-S. 2018 . Spatial variation in habitat quality for juvenile flatfish in the southeastern Bering Sea and its implications for productivity in a warming ecosystem . Journal of Sea Research , 139 : 62 – 72 . Google Scholar Crossref Search ADS WorldCat Published by International Council for the Exploration of the Sea 2019. This work is written by US Government employees and is in the public domain in the US. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Contrasting the variability in spatial distribution of two juvenile flatfishes in relation to thermal stanzas in the eastern Bering Sea JO - ICES Journal of Marine Science DO - 10.1093/icesjms/fsz180 DA - 2020-05-01 UR - https://www.deepdyve.com/lp/oxford-university-press/contrasting-the-variability-in-spatial-distribution-of-two-juvenile-ODTbg1m1AI SP - 1 VL - Advance Article IS - DP - DeepDyve ER -