key: cord-0041662-docw2lwy
authors: Busby, M. S.; Duffy‐Anderson, J. T.; Mier, K. L.; De Forest, L. G.
title: Spatial and temporal patterns in summer ichthyoplankton assemblages on the eastern Bering Sea shelf 1996–2007
date: 2014-03-11
journal: Fish Oceanogr
DOI: 10.1111/fog.12062
sha: 98140f37198777fcf82b0240735c73cd7f5571d6
doc_id: 41662
cord_uid: docw2lwy

Larval and early juvenile fishes were sampled from the eastern Bering Sea (EBS) shelf from 2001 to 2005, and in 2007. Data from these collections were used to examine spatial and temporal patterns in species assemblage structure and abundance. The years 2001–2005 were unusual because the EBS water temperature was ‘warm’ compared with the long‐term mean temperature. In contrast, 2007 was a ‘cold’ year. The abundance of the five most numerous taxa at 12 stations common to all years sampled (1996–2005, 2007) were significantly different among years. Larval and early juvenile stage Theragra chalcogramma (walleye pollock), a commercially important gadid, were by far the most abundant fish in all years. Bottom depth alone best explained assemblage structure in most years, but in others, bottom depth and water column temperature combined and percent sea‐ice coverage were most important. Abundance of T. chalcogramma larvae increases with water column temperature until 5°C and then becomes level. Higher abundances of Gadus macrocephalus (Pacific cod) larvae occur in years with the greatest percent sea‐ice cover as indicated by GAM analysis. Larvae of Lepidopsetta polyxystra (northern rock sole) increase in abundance with increasing maximum wind speed, but decrease at a later date during the last winter storm. The data are consistent with the hypothesis that oceanographic conditions, specifically water temperature and sea‐ice coverage, affect the spatial and temporal pattern of larval abundances. In general, ichthyoplankton species assemblages can be important early indicators of environmental change in the Bering Sea and potentially other subarctic seas as well.

The distribution of planktonic organisms in marine environments has been central to numerous studies throughout oceanic and coastal areas of the world or 'large marine ecosystems' (Sherman, 2000; Sherman et al., 2007) . Planktonic fish eggs and larvae, commonly referred to as ichthyoplankton, have been collected in regional surveys for many purposes including studies of ecosystem functions, stock assessment, and investigations of recruitment processes (Kendall and Duker, 1998; Kendall, 2000) .

It is well known that the complexity of species assemblages and biodiversity is related to ecosystem health and functioning (Magurran, 2004) , and that links between ecosystem variations and shifts in species richness have been shown in marine systems across the globe (Meyers and Worm, 2003; Micheli and Halpern, 2005; Wilson et al., 2006; Cheung et al., 2009) . Many of these studies have focused on the interactions among adult members of the community, but the study of species associations among earlier life stages, such as larvae and juveniles, may also be important harbingers of phase shifts in ecosystem dynamics. In particular, ichthyoplankton are greatly influenced by dynamic physical processes, including oceanographic features such as bathymetry (Doyle et al., 2002) , eddies (Atwood et al., 2010) , currents and fronts (Norcross et al., 2010; Siddon et al., 2011) , as well as meteorological features such as precipitation and freshwater discharge (Boeing and Duffy-Anderson, 2008) , wind variation (Rodriguez et al., 2011) , and temperature (Genner et al., 2010) , which can lead to shifts in community structure, changes in biotic control, and variable species interactions. The sensitivity of early life stages to these and other forcing phenomena potentially make marine fish larvae important sentinels of regional climate variation and ecosystem change.

The Bering Sea is one of the most productive ecosystems in the world and accounts for approximately 5% of the world fish and shellfish catch and over 40% of all US finfish and shellfish landings (Stabeno et al., 2012a) . The eastern Bering Sea (EBS) has been shown to respond rapidly to climate (Napp and Hunt, 2001) . Increasingly, and in light of recent climate variability, there is concern over the present status of the Bering Sea and the future of its ecosystem health (Grebmeier et al., 2006) . Studies that examine processes that might be driving ecosystem regime shifts (Niebauer and Day, 1989; Stabeno et al., 2001 Stabeno et al., , 2007 , what the environmental effects are on the biota (Brodeur et al., 1999; Kotwicki et al., 2005) , what the consequences are for biodiversity and species interactions, and how long effects may be expected to last (Benson and Trites, 2002; Duffy-Anderson et al., 2006) are critical to the evaluation of potential impacts of climate change.

The extent and concentration of sea-ice coverage and timing of retreat directly influence the dynamics of primary and secondary production, which in turn have substantial effects on energy flow and ecosystem function in the EBS . Winter sea-ice concentration in the EBS has declined since the mid-1970s, but some marked fluctuations in yearly mean water column and sea surface temperatures (SST), manifested as alternating periods of variable duration, have been noted and classified as either 'warm' or 'cold' years (Stabeno et al., 2012a ). Depthaveraged spring and summer temperatures are closely related to the previous winter's temperatures and ice cover. Variations in depth-averaged temperatures and timing of sea-ice retreat directly influence the timing and location of the spring phytoplankton bloom and the availability, abundance, and species composition of zooplankton, particularly copepods, which are the primary components in the diets of fish larvae and juveniles (Hunt et al., 2011) . During the previous decade, an unusual event occurred, characterized by five consecutive 'warm' years (2001-2005) followed by a 'normal' or transition year (2006) and four consecutive 'cold' years (2007-2010) (Stabeno et al., 2012a) . This provides a unique opportunity to study the biological and physical processes of the EBS ecosystem in two different climatic 'regimes.' The temperatures recorded during this recent warm period are thought to be characteristic of those predicted by 2050 under a climate warming scenario (Stabeno et al., 2012b) .

The influence of ice retreat timing and related consequences to ecosystem dynamics and function persist into the summer. Species composition, relative abundance, and length of larval and juvenile fishes in Bering Sea ichthyoplankton assemblages during summer are of particular interest as this is a period of ontogenetic transition for many species. For example, larvae of winter and spring spawning flatfish species are undergoing transformation and eye migration is nearly complete. Age-0 Theragra chalcogramma (walleye pollock), a commercially important gadid fish, form scales and begin to resemble adults in morphology and behavior (Brown et al., 2001) , which includes forming schools and feeding at night (Merati and Brodeur, 1996; Sogard and Olla, 1996; Olla, 1997, 1998) . It is during this critical stage of development that many demersal fish species make the transition from pelagic to benthic habitats and are particularly vulnerable to mortality from predation and other causes (Toole et al., 1993) . Survival of age-0 T. chalcogramma is linked to biochemical condition prior to onset of their first winter Heintz et al., 2013) . Changes in the timing of transitions between early life history stages due to substantial shifts in prevailing environmental conditions may have significant consequences for survival and recruitment of certain species and the structure of species assemblages.

Our objectives in this study were to (i) describe midwater assemblages of larval and juvenile fishes in the EBS in summer for 2001-2005, and 2007; (ii) relate observations of fish abundances and lengths to oceanographic conditions; (iii) examine the interannual stability of species assemblages using data from the present study (2001) (2002) (2003) (2004) (2005) 2007) and data from a prior investigation Duffy-Anderson et al., 2006) , which together span a series of warm and cold regimes over the EBS shelf, and (iv) compare results of our analyses of the 2001-2005 and 2007 ichthyoplankton data with those from 1996-2000. To accomplish these objectives we adopted a two-step approach. First, we undertook a within-year examination of species assemblages during the most recent years of collection (2001) (2002) (2003) (2004) (2005) 2007) and made quantitative and qualitative comparisons to assemblage structure from 1996-2000. Next, we conducted a time-series analysis over the entire data set (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) to examine interannual variations in species lengths and abundances.

Ichthyoplankton samples used in this study were collected from the EBS shelf during a series of cooperative research cruises between scientists at the Alaska Fisheries Science Center (AFSC), Seattle, WA, USA, and the Graduate School of Fisheries, Hokkaido University, Japan, on the T/S Oshoro Maru during the years 2001-2005 and in 2007 . In July (and sometimes August) of each year (Table 1) , a grid of stations on the EBS shelf was occupied, extending mostly eastward and south from the Pribilof Islands (Fig. 1) . Conductivity-temperature-depth (CTD) casts were made immediately prior to fish sampling at nearly every station in each year. A modified beam trawl (MBT), with a 5-m 2 mouth opening fitted with a 3 9 2 mm oval mesh net and a 1-mm mesh cod end, was towed obliquely from 200 m (or 10 m off bottom if shallower) to the surface. Net depth was measured using an acoustic netsonde. The MBT was designed to collect late larval and early juvenile stage fishes in the midwater and is similar to the Methot trawl (Methot, 1986) . The MBT differs from the Methot trawl in not having a depressor and weight is contained within the tubular frame. A flowmeter was fitted into the mouth of the net frame to determine volume filtered. Samples were preserved in 5% formaldehyde-seawater solution buffered with sodium borate. Because of the large mesh size of the MBT, preflexion and early flexion stage larvae of all species were likely not quantitatively sampled. Supplemental surface and bottom temperature data were obtained from AFSC annual bottom trawl surveys of the EBS shelf in summer (June-August). Only temperature data collected from within the same geographic extent as the present study were used. Field collection methods generally follow those of Duffy-Anderson et al. (2006) .

Preserved samples were sorted and fishes identified to the lowest taxonomic level possible at the Plankton Sorting and Identification Center in Szczecin, Poland. Taxonomic identifications were verified at the AFSC. Some fish were categorized as taxonomic groups (e.g., Atheresthes spp., Sebastes spp.) due to limitations associated with identifying larval stages to the species level. Fish were measured for standard length (SL) to the nearest 1.0 mm. The separation point between the larval and juvenile stages for Theragra chalcogramma and Gadus macrocephalus is 25.0 mm SL based on the size at the transformation of T. chalcogramma determined by Brown et al. (2001) . Catch in each tow was converted to catch per 10 m 2 of sea surface area.

Two data sets were used. The first from field sampling conducted in 2001-2005 and in 2007 was analyzed to describe species associations and to relate patterns to predominant oceanographic conditions. The second data set combined previously published data from (Duffy-Anderson et al., 2006 with the 2001-2005 and 2007 data to derive a time series (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007 ) that could be used to examine interannual variability in assemblage structure, species abundance and length, and environmental associations. Environmental data Bottom depth, temperature, and salinity were recorded at each sampling station using a CTD data recorder. Sea ice data were collected using the special sensor microwave interferometer (SSMI) and obtained from the National Snow and Ice Data Center (NSIDC, http://nsidc.org/). Data were extracted into nominally diagonal 25 9 25 km regions and reported as percent ice coverage for each. Definitions of 'warm' and 'cold' years are based on depth averaged temperature criteria reported by Stabeno et al. (2012a) . The wind velocity data were derived using algorithms from the NCEP/NCAR Reanalysis data set (Kalnay et al., 1996) . This data set was produced using a numerical weather prediction model with the assimilation of available observations. Ladd and Bond (2002) have shown that it provides an accurate characterization of the winds in the Bering Sea. The present project incorporates daily averaged surface wind speeds for the years 1996-2007 at the gridpoint locations 56.2°N, 165.0°W and 56.2°N, 168.8°W. The data were used to determine the timing (date) of the first storm event after ice retreat (last winter storm) similar to that shown by . Criteria for defining a 'storm' includes a minimum 3-day period of winds exceeding 5 ms À1 with at least one day being >10 ms À1 . The date assigned to each storm was the day of the year (DOY) at the midpoint of the event. The maximum wind velocity reported from each storm event was used as a variable to represent the force of mixing in the upper water column.

Spatial analyses (2001) (2002) (2003) (2004) (2005) 2007) Methods used to evaluate species composition for 1996-2000 data (Duffy-Anderson et al., 2006) were applied to the 2001-2005 and 2007 data. Only species occurring in at least 4% of the samples for a particular year were included. Removal of rare species was deemed appropriate because these were not likely to contribute significantly to the broader spatial and temporal patterns. Separate hierarchical clustering analyses for species and stations (using abundance as catch per 10 m 2 ) were performed for each year using the Bray-Curtis dissimilarity coefficient as the distance measure. A hierarchical agglomerative cluster analysis was used, specifically the flexible beta clustering method with beta = À0.25. A fourth-root transformation was applied to the abundance data before calculating the Bray-Curtis dissimilarity coefficient for clustering stations (samples) to more equally weight the contribution of rare and abundant species.

Interannual analyses (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) Differences in species composition across all years (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) were examined using a two-way analysis of similarity (ANOSIM), using year and station as factors. As with the cluster analysis, the matrix of Bray-Curtis dissimilarity coefficients calculated on fourth-root transformed catch (10 m À2 ) was used as input. ANOSIM is a nonparametric, multivariate analysis that uses permutation techniques to test for differences in species composition among years, analogous to a univariate analysis of variance (ANOVA). When the overall test was found to be significant, permutational pairwise multiple comparison tests were conducted. Wherever there were significant differences (P ≤ 0.05), a similarity percentages (SIMPER) analysis (Clarke and Warwick, 2001 ) was used to calculate each species contribution to the Bray-Curtis dissimilarity, therefore defining the discriminating species that best explain each significant difference.

Variations in abundances and lengths among all years were investigated for the five most abundant taxa. Since there were slight variations in the geographic area sampled in each year a subset of 12 stations spanning the middle and outer shelf were identified as having been sampled in all study years, and variations in abundance and length among the years were evaluated from this 'core' group of 12 stations (Fig. 1 ). These analyses were conducted for each of the five most common taxa within this core group of stations; T. chalcogramma, Hippoglossoides elassodon, G. macrocephalus, Lepidopsetta polyxystra, and Atheresthes spp. A mixed model analysis of length data, weighted by the haul catch per unit effort (CPUE) was used to test for length differences among years, and ANOVAs were used to test for differences in abundance, using station location as a block. These were followed by Tukey multiple comparison tests to see which years were significantly different.

Environmental analyses (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) A quantitative biotic-environmental (BIO-ENV) analysis was applied (Clarke and Warwick, 2001) to each year (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) to determine which physical variables were most related to the species composition based on abundance. The physical variables were bottom depth, integrated water column temperature, bottom temperature, salinity, and percentage of area covered by ice in winter. To investigate the relationships of environmental variables with interannual differences in abundance of the five most abundant taxa, a generalized additive model (GAM) was employed based on the 'core' stations ( Fig. 1) .

The variables used were water column temperature, salinity, percent sea-ice coverage, date (DOY) of last winter storm, and maximum wind speed during last winter storm. The GAM was used instead of multiple linear regression to allow for possible nonlinear relationships as well as to relax the normality assumption. The best fit GAM was determined by checking all possible combinations of the environmental variables and selecting the one with the lowest GCV (general crossvalidation) as well as the one containing variables with P-values <0.05. composition (2001-2005, 2007) A total of 18 290 larval and juvenile fishes comprising at least 43 taxa representing 17 families was collected during the study period. These fish were predominantly in the postflexion and transformation stage of development. Taxonomic richness was lowest in 2002 and 2003 with only 17 and 18 taxa collected, respectively, compared with 21-28 in other years ( Table 2 ). The family Pleuronectidae (flatfishes) was represented by the greatest number of taxa (10), followed by Cottidae and Agonidae (7). Taxonomic richness within certain families was variable throughout the study period. Particularly noteworthy is the family Cottidae where only one of the seven taxa reported for the 6 yr studied were caught in 2002 (Hemilepidotus hemilepidotus) and in 2004 (Icelus spp.). Richness of pleuronectids was lowest in 2003, when only three taxa (Atheresthes spp., H. elassodon, and L. polyxystra) were collected.

The gadids T. chalcogramma and G. macrocephalus were caught every year, as were Liparis spp., Lumpenus maculatus, Ammodytes hexapterus, Atheresthes spp., H. elassodon, and L. polyxystra. Seven taxa were collected for 5 of the 6 yr; Sebastes spp., H. hemilepidotus, Bathyagonus alascanus, Podothecus veternus, Zaprora silenus, L. bilineata, and Reinhardtius hippoglossoides. (Table 3) . Hippoglossoides elassodon, G. macrocephalus, L. polyxystra, and Atheresthes spp. made up the remainder of the five most common taxa in the order listed, but when combined only accounted for 5.6% of the total catch. Spatial patterns (2001) (2002) (2003) (2004) (2005) 2007) Results In 2001, the first of five consecutive 'warm' years, the assemblage structure consisted of five species groups and three station groups (Fig. 2) . A relatively cohesive outer domain (100-200 m) (station group 1; Fig. 3 ) assemblage consisted of moderate abundances of T. chalcogramma (species group D) with lower abundances of C. aleutensis, G. macrocephalus, H. elassodon, and Z. silenus (species group E). The middle domain (50-100 m) assemblage (station group 3) was dominated by T. chalcogramma with L. polyxystra. For 2001-2004 station, group 2 straddled the middle and outer domain and although T. chalcogramma was also in high abundance as in station group 3, there were no unique species (not present in the other years) present or distinct patterns (Table 4, Fig. 3 ).

In 2005, the last of five consecutive 'warm' years, the highest number of taxa (n = 28, Table 2 ) were present. Cluster analysis revealed eight species groups, but only two were well defined by having one of the two dominant taxaspecies groups C and E (Table 4 , Fig. 4) . Species group C consisted of G. macrocephalus in moderate abundance with Z. silenus, Icelinus borealis, and Liparis spp. also present. This group was mostly associated with station group 3 in the outer domain and near the Pribilof Islands (Fig. 3) . Station group 4 (Figs 3 and 4 ) over the middle domain had the highest abundances of T. chalcogramma. The other Table 2 . Fish taxa present in at least 4% of the samples for each year. 'X' indicates whether that taxon was collected in a particular year. The total number of taxa collected in each year is given at the bottom (total includes rare taxa that were removed before the analysis). Flathead Sole X X X X X X taxa present were in low abundance and species in groups A, B, D, and F were absent. Station group 3, mostly in the outer domain and Pribilof Island region, had greater richness and all taxa in species group E were present at all stations with the exceptions of H. elassodon and T. chalcogramma, which were each absent at one station. Station group 1 and an outlier station located off the edge of the continental shelf were represented by the deepwater taxa Coryphaenoides spp., L. schmidti, S. leucopsaurus, and Sebastes spp., and G. macrocephalus, M. zonurus, and B. nigripinnis, which are not typically considered deep water species, were also present.

In 2007 there were four species and three station groups (Table 4 , Fig. 5 ) with the clearest distinction of species and station groups reflecting both bathymetry and temperature patterns (Figs 3 and 5) . It was also the first truly 'cold' year since 1999 and had the lowest abundance of T. chalcogramma over the entire 1996-2007 time series (Fig. 6 ). There were no stations off the continental shelf break as in previous years and thus no deepwater species were collected. Station group 1 located in the outer domain comprised only taxa from group D that included Atheresthes spp., Dasycottus setiger, G. macrocephalus, and H. hemilipidotus, in addition to T. chalcogramma, which was caught at every station. Station group 2 in the middle domain comprised all seven taxa in species group C, of which all but L. bilineata were unique to station group 2. Species group B, composed of H. elassodon and R. hippoglossoides, was also unique to station group 2. The three taxa in species group A and A. hexapterus were only present at the two stations in station group 3 straddling the boundary of the inner and middle domains. Three species, U. olriki, R. jordani, and L. proboscidea, were only collected in 2007 (Table 2) .

In general, with the exception of 2003 and 2007, there was a deepwater assemblage present, usually consisting of the bathylagids B. pacificus or L. schmidti, macrourids Coryphaenoides spp., scorpaenids Sebastes spp., or pleuronectids Atheresthes spp. Another noteworthy result from the cluster analyses was the very different taxa each year in the species groups that included T. chalcogramma (Table 4 ). In 2001, for example, L. polyxystra was the second most frequently occurring and abundant species (and thus most closely associated with T. chalcogramma: Fig. 2) , whereas H. elassodon took this position in 2002 and 2005 (closely followed by Atheresthes spp.). In 2003, T. chalcogramma dominated species group B, occurring at all but two stations; the other species in this group, A. hexapterus and L. maculatus, occurred at only two and one stations respectively. In 2004, T. chalcogramma was closely affiliated with G. macrocephalus, which although in low abundance, co-occurred at several stations. Dasycottus setiger was the species most closely associated with T. chalcogramma in 2007.

Interannual comparisons (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) The ANOVA and multiple comparison test results showing differences in abundance among all 11 yr of sampling at the 'core' stations indicated that T. chalcogramma was most abundant in 2001 but with no significantly higher abundance than in 2000 and 2004, and lowest in abundance in 1998, 2002 and 2007 (Table 5 , Fig. 6 ). Hippoglossoides elassodon was present in high abundance in 1998, 2001, 2002, 2003, 2004, and 2005 and in low abundance in , 1997 , 1999 , 2004 . Gadus macrocephalus abundance was highest in 1996 , 2003 , 2004 . Fig. 7) . Although the length of T. chalcogramma in 2004 was significantly less than in 1996 and 2005, the length was significantly larger than all other years except 1998 (P < 0.003). Significant differences in length were also detected among years for G. macrocephalus, Atheresthes spp., H. elassodon, and L. polyxystra, but lacked any consistency.

Environmental patterns (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) The BIO-ENV procedure conducted within each year (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) 2007) determined which environmental Table 4 . Species groupings as identified by cluster analyses. Proximity of letter groups signifies distance between species clusters in dendograms for each year (refer to Fig. 2, 2001 for comparison) . Order of taxa present within letter groups in each year indicates proximity of adjacent taxa within a cluster. Dash indicates cluster group not present. Refer to Table 2 factors were most consistent with the ichthyoplankton assemblage structure. Bottom depth was usually, but not always, the environmental factor most correlated with taxonomic composition (Table 6 ). In , 1997 In -2003 , bottom depth alone was most correlated with assemblage structure, but in 1998, 1999, and 2007 bottom depth and water column temperature were equally correlated. In 2004, bottom depth and salinity were the environmental variables most highly correlated with assemblage structure, and in 2000, percent sea-ice coverage was correlated with ichthyoplankton assemblage structure. Date (DOY) and strength of last winter storm were not correlated with assemblage structure. A correlation analysis of environmental variables indicated that average water column temperature and bottom temperature were correlated with one another; therefore, bottom temperature was not included in the GAM analyses. GAM analyses conducted to evaluate potential influences of environmental factors on abundances of the five most common taxa at the 12 'core' stations over all years sampled, indicated that abundance of T. chalcogramma larvae increased with increased mean water column temperature until about 5°C, at which point it leveled off (P = 0.003; Fig. 8a ). Percent sea-ice cover was significant for G. macrocephalus (P = 0.02), with the abundance of larvae increasing with increased percent sea-ice cover (Fig. 8b ). Maximum wind speed had a significant positive linear effect for L. polyxystra (P = 0.005) when combined with the date (DOY) of the last winter storm (abundance of larvae increasing with wind speed, Fig. 9a ), whereas abundance decreased with later date of the last winter storm (Fig. 9b) .

We hypothesized that ichthyoplankton species associations would be sensitive to environmental perturbations, and anticipated that we would observe differences in assemblage structure related to anomalously warm or cold thermal conditions over the shelf. Indeed, we observed changes in ichthyoplankton assemblage structure in cold years (1997, 1999, 2007) relative to warm (1998 ) and 'average' years (1996 ) (Duffy-Anderson et al., 2006 present study) . A contributing factor to these changes were shifts in species composition consistent with ingress of warm and cold water species in warm and cold years.

In comparing the 10 most abundant taxa from the 1996-2000 sampling (alternating warm and cold years) with 2001-2005 and 2007 (five consecutive warm years and one cold year), we observed replacement of three taxa ( Table 7) . Two of these, Sebastes spp. and R. hippoglossoides, are deepwater or slope taxa, the latter being a pleuronectid with a circumpolar distribution that might be considered an indicator species affiliated with colder Arctic water masses (Mecklenberg et al., 2011) . The abundance of R. hippoglossoides larvae declined somewhat during the warm period, suggesting a shift in prevailing oceanographic conditions (Fig. 10a) . More frequent occurrences during the warm years (2001) (2002) (2003) (2004) (2005) of larval L. bilineata, a pleuronectid with a more southerly range, may be indicative of a northern range expansion under conditions of rising temperatures (Fig. 10b) . It should also be noted, however, that occurrences of adult L. bilineata in AFSC groundfish surveys are limited to only five stations immediately north and east of Unimak Island (Fig. 1 ) and the individuals collected were relatively small and probably not of spawning age (D. Stevenson, AFSC pers. com.). However, it is possible that larvae originating from adult L. bilineata spawning in the Gulf of Alaska are transported northward through Unimak Pass and onto the southeast Bering Sea shelf, as the prevailing currents are often very strong here with considerable seasonal and interannual variability Duffy-Anderson et al., 2006) . It is not surprising that our 10 most abundant taxa apparently respond differently in terms of larval/juvenile abundance (see Table 7 ) to environmental influences, as this is most likely related to the evolutionary history of each taxon.

Moreover, we observed that oceanographic conditions, notably temperature and sea-ice, have a significant effect on the larval abundances of two Figure 6 . Interannual variations in abundance AE standard error (SE) of the five most abundant taxa collected from the 'core' stations. Note differences in abundance scales between species (number 10 m 2 ). -Anderson et al., 2006) was that assemblages in 1998 were significantly different from other years examined, and that a short-term fluctuation in oceanographic conditions, the unusually strong 1997-1998 El Niño, appeared to disrupt predominant species assemblage patterns over the continental shelf, although effects were transitory. Lower larval abundances were also observed during the El Niño, with a return to previous levels the following year. Indeed, unique patterns were observed that appeared to be related to thermal conditions. Revisiting the 1996-2000 data examined in the initial study provided the opportunity to update some (Cartwright, 2009 ). This reduced the number of occurrences of Icelus spp. below the 4% threshold required for inclusion in the analysis. Psychrolutes spp. was included in the 5% group in the initial study, but only occurs in 2.6% of the 1996-2000 hauls in our taxonomically updated data set. In addition, Podothecus veternus was one of the agonid species collected in the recent sampling, not P. acipenserinus as previously reported. Although these specimens previously identified in the initial study were not re-examined, it is likely that they are mostly, if not all, P. veternus.

Sampling in additional 'cold' years after 2007 would have allowed us to gain a more in-depth understanding of assemblage structure and variability under these conditions. It appears from our observations in 2007 that assemblage structure is more cohesive and boundaries between species and station clusters (groups) and oceanographic (thermal) fronts are more distinct than in 'warm' years. During the five consecutive 'warm' years (2001) (2002) (2003) (2004) (2005) , the early ice retreat likely led to low standing stocks of large zooplankton species, which may have limited overwinter survival of age-0 T. chalcogramma by depriving them of important energy-rich prey resources in summer and autumn (Hunt et al., 2011) . This hypothesis is based on more recent observations determining that cold conditions lead to greater production and availability of large higher-energy-content zooplankton such as the euphausiids Thysanoessa raschii and T. inermis and the copepod Calanus marshallae . Although 2007 was the first 'cold' year after a transition (2006, an average year, unsampled) it had the lowest age-0 T. chalcogramma abundance over the entire time series. In addition, overall zooplankton availability in 2007 was high, but survival of the age-0 T. chalcogramma was low compared with subsequent cold years 2008-2010 . One possible explanation for this is that the standing stock of larger zooplankton species with high lipid content had not yet fully recovered after several years of poor production to replace the smaller lipid poor species that were very abundant. Overwinter survival is also important for recruitment of these zooplankton species in the following year.

Unfortunately, this study could not determine the actual source of interannual variations in larval Rank 1996 Rank -2000 Rank 2001 Rank -2005 Rank , 2007 1 Theragra chalcogramma Theragra chalcogramma 2

Hippoglossoides elassodon Hippoglossoides elassodon 3

Gadus macrocephalus Gadus macrocephalus 4

Atheresthes spp. Lepidopsetta polyxystra 5

Lepidopsetta polyxystra Ammodytes hexapterus 6

Ammodytes hexapterus Atheresthes spp. 7

Lumpenus maculatus Liparis spp. 8

Reinhardtius hippoglossoides Podothecus veternus 9

Sebastes spp.

Zaprora silenus Lumpenus maculatus abundances of the taxa studied. Differences could be due to variations in larval supply, mortality, advection, and/or availability of nursery habitat, among other things. Moreover, observed trends in larval and juvenile abundance do not necessarily correlate with trends in age-1 recruitment estimates (Ianelli et al., 2010) , suggesting that other factors also play a role in recruitment. Nevertheless, our results show that the early life history stages of marine fishes do respond to environmental perturbations relatively quickly, in terms of both changes in their abundances and/or changes in their species associations. The data we presented here encompass late larval, transformation, and early juvenile stages, but the influences of gradual environmental change or irregular 'disruptions' on spawning behavior (timing), hatching, and survival of earlier stages must also be considered as there are remnant or linked effects to later stages. Perhaps the most important conclusion from this study is that an 11-yr time series is not long enough to fully understand the underlying mechanisms that influence change in ichthyoplankton assemblages and adult fish populations. As such, ichthyoplankton metrics remain useful harbingers of environmental variability and ecosystem change.

Influence of mesoscale eddies on ichthyoplankton assemblages in the Gulf of

Ecological effects of regime shifts in the Bering Sea and eastern North Pacific Ocean

Ichthyoplankton dynamics and biodiversity in the Gulf of Alaska: responses to environmental change

Evidence for a substantial increase in gelatinous zooplankton in the Bering Sea, with possible links to climate change

Transformation from the larval to juvenile stage in walleye pollock (Theragra chalcogramma) -the relationship between otolith and osteological development

Description of early life history stages of the northern sculpin (Icelinus borealis Gilbert) (Teleostei: Cottidae)

Projecting global marine biodiversity impacts under climate change scenarios

Change in Marine Communities: An Approach to Statistical Analysis and Interpretation

Climate change in the southeastern Bering Sea: impacts on pollock stocks and implications for the oscillating control hypothesis

Regional variation in springtime ichthyoplankton assemblages in the northeast Pacific Ocean

Spatial and temporal patterns in summer ichthyoplankton assemblages on the eastern Bering Sea shelf

Temperature-driven phonological changes within a marine larval fish assemblage

A major ecosystem shift in the northern Bering Sea

Climate related changes in the nutritional condition of young-of-the year walleye pollock (Theragra chalcogramma) from the eastern Bering Sea

Climate change and the control of energy flow in the southeastern Bering Sea Ecosystem

Climate change and control of the southeastern Bering Sea pelagic ecosystem

Climate impacts on eastern Bering Sea foodwebs: a synthesis of new data and an assessment of the Oscillating Control Hypothesis

Eastern Bering Sea Walleye Pollock. pp 53-156 In: 2010 North Pacific Groundfish Stock Assessment and Fishery Evaluation Reports for 2011 North Pacific Fishery Management Council, 605 W

The NCEP/ NCAR 40-year reanalysis project

Status of recruitment studies of northeast Pacific fishes

The development of fisheries oceanography in the United States

Variation in the distribution of walleye pollock (Theragra chalcogramma) with temperature and implications for seasonal migration

Evaluation of the NCEP-NCAR reanalysis in the northeast Pacific and the Bering Sea

Measuring Biological Diversity

Biodiversity of arctic marine fishes: taxonomy and zoogeography

Feeding habits and daily ration of juvenile walleye pollock, Theragra chalcogramma, in the western Gulf of Alaska

Frame trawl for sampling juvenile fish. CalCOFI Rep

Rapid worldwide depletion of predatory fish communities

Low functional redundancy in coastal marine assemblages

Anomalous conditions in the south-eastern Bering Sea 1997: linkages among climate, weather, ocean, and biology

Causes of interannual variability in the sea ice cover of the eastern Bering Sea

Demersal and larval fish assemblages in the Chukchi Sea

Composition and daytime vertical distribution of the ichthyoplankton assemblages in the Central Cantabrian Sea shelf, during summer: an Eulerian study

Altered search speed and growth: social versus independent foraging in two pelagic juvenile fishes

Effect of light on juvenile walleye pollock shoaling and their interaction with predators

Marine ecosystem management of the Baltic and other regions

Global Applications of the Large Marine Ecosystem Concept

Community-level response of fish larvae to environmental variability in the southeastern Bering Sea

Conceptual model of energy allocation in walleye pollock (Theragra chalcogramma) from age-0 to age-1 in the southeastern Bering Sea

Diel patterns of behavior in juvenile walleye pollock Theragra chalcogramma

On the temporal variability of the physical environment over the south-eastern Bering Sea

Transport through Unimak Pass

On the recent warming of the southeastern Bering Sea shelf

Comparison of warm and cold years on the southeastern Bering Sea shelf. Deep-Sea Res

A comparison of the northern and southern shelves of the eastern Bering Sea and some implications for the ecosystem. Deep-Sea Res

Relationships between otolith microstructure, microchemistry, and early life history events in Dover Sole, Microstomus pacificus

Multiple disturbances and the global degredation of coral reefs: are reef fisheries at risk or resilient?

The authors thank the officers, crew, and students (cadets) of the T/S Oshoro Maru. Sigrid Salo and Nick Bond assisted by obtaining sea-ice coverage and wind data. Miriam Doyle, Ann Matarese, Jeff Napp, and Phyllis Stabeno provided guidance and helpful suggestions throughout the study and on earlier drafts of the manuscript. Deborah Blood provided editorial assistance. The findings and conclusions in this paper are those of the authors and do not necessarily represent the views of the National Marine Fisheries Service, NOAA. Reference to trade names does not imply endorsement by the National Marine Fisheries Service, NOAA.