Introduction

Polar cod (Boreogadus saida), also known as Arctic cod, is considered a key species in Arctic ecosystems (Andriashev 1954; Bradstreet et al. 1986; Hop and Gjøsæter 2013). It has a circumpolar distribution and is normally found in the eastern and northern Barents Sea, where it plays a key role as forage fish (Gjøsæter 2009; Ajiad et al. 2011). The stock structure of polar cod in the Barents Sea and adjacent areas is not well known, but based on new modelling work (Huserbråten et al. 2019), and past anectodal evidence (Ponomarenko 1968) the stock appears to be split in at least two components, one spawning in the southeastern parts of the Barents Sea known as Pechora Sea and one spawning at annually variable latitudes east of Svalbard. The current knowledge on the geographical distribution of polar cod in the Barents Sea is based on surveys that primarily targeted capelin, and areas north and east of the usual distribution of capelin have thus not been covered. The population may, therefore, extend beyond the survey area, particularly along the northern and eastern edges (Ajiad et al. 2011). In recent years, a reduction of stock size from 1 to 2 million tonnes to less than half of that, weak recruitment and a displacement towards the northern and eastern parts of the Barents Sea have been observed, which are often attributed to climatic changes and a “borealization” of the Barents Sea (Fossheim et al. 2015; Eriksen et al. 2017). However, as pointed out by Gjøsæter (1995), the survey-derived stock estimates should be treated as a relative index of abundance, as the survey was originally developed to assess capelin abundance and does not fully encompass the area occupied by polar cod. The distribution of polar cod continues eastwards into the Kara Sea and further along the Siberian Shelf (Ponomarenko 1968), and this may in fact represent a continuation of the distribution area of the stock spawning in the Barents Sea.

The annual surveys, in which polar cod was covered, have been conducted in autumn (mostly September) since 1972, but have undergone major changes over the years. In the early period, the survey was a capelin and polar cod survey, then a pelagic fish survey, then a multispecies survey, and finally, from 2004, an “ecosystem survey” known as the Barents Sea Ecosystem Survey (BESS) (Gjøsæter 2011; Eriksen et al. 2018).

Sampling protocol and survey description can be found in Toresen et al. (1998), Michalsen et al. (2013) and Eriksen et al. (2018). For the most part, data over the total distribution area of polar cod were aggregated to form a stock size estimate for the total stock. In addition to the acoustic abundance estimates of polar cod age 1 and older, an abundance index of age 0 polar cod is constructed from dedicated trawl hauls with a pelagic trawl towed in the upper 60 m (Eriksen et al. 2009). Before 2004, this index was based on a 0-group survey conducted in late August, but in 2004 the 0-group investigations became part of the BESS (Eriksen et al. 2018). The timing of the 0-group investigations did not change much.

In the wake of the increased attention, climate change has attained in recent years (IPCC 2014, 2019), and the increased sea temperatures that have already been observed in the Barents Sea and other arctic areas (Smedsrud et al. 2013; ICES 2017, 2018), the fate of polar cod and other true arctic species has become a theme of several research projects and research articles (Christiansen 2017; Kunz et al. 2018; Bender et al. 2018; Koenker et al. 2018; Laurel et al. 2018). The motivation for studying polar cod in the climate change context is driven by the fact that this species has a circumpolar distribution and in many of the areas where it occurs, it is the dominant pelagic species (Ajiad et al. 2011; Bouchard et al. 2017), playing a dominant role in those ecosystems (Lowry and Frost 1981; Bradstreet et al. 1986; Orlova et al. 2009; Hop and Gjøsæter 2013; Mueter et al. 2016).

The polar cod stock in the Barents Sea is probably among the largest occurrences of polar cod, and has been studied for a century or more, first and foremost by Russian scientists, but also since 1960s by Norwegians. The Barents Sea is already experiencing increased sea temperatures and, given that it is predicted to experience faster and more profound climatic changes in the future than most other areas (Serreze and Barry 2011; Smedsrud et al. 2013), studies that investigate how the life history and ecology of polar cod in the Barents Sea are impacted by climate change may provide insight into what may happen in other arctic areas. In addition to climate change, anthropogenic activities like petroleum and mineral exploration, transportation and tourism that may follow from reduced sea ice and increased access to hitherto inaccessible areas has increased the interest in studying arctic biota inhabiting the northern Barents Sea.

The goal of this research paper is (1) to use the output of a biophysical individual-based model (IBM) of polar cod to characterize the environmental and developmental properties (spawn time, hatch time, ice cover, temperature) of the early life history of individuals that reached the 0-group stage at the time and place of observations and test if they varied between eastern and western spawning locations, and (2) to use a statistical linear regression model to examine if and how ice concentration, ice breakup time, maximum temperature, and spawning stock biomass relate to modelled larval recruitment success to the 0-group, hereafter defined as spatial and temporal match between simulated individuals from potential spawning grounds and the observed 0-group stage. Also, when discussing the recruitment success of larvae, depending on context, this specifically means the successful recruitment of eggs and larvae to the 0-group, as distinct from the traditional definition as recruitment to the adult stock.

Increased insight into how reduction of sea ice will impact early life stages of polar cod is needed to better understand how the ongoing climate change may affect recruitment and thereby the whole population dynamics of this ecologically important fish species.

Short description of study area

The Barents Sea (Fig. 1) is a high-latitude, shallow (average depth 220 m) shelf sea bounded in the north and west by the continental shelf edge towards the Arctic Ocean (at about 80–81°N) and the Norwegian Sea (at about 10–15°E). The hydrographic situation is mainly affected by three water masses; Atlantic, Coastal and Arctic waters (Loeng 1991; Smedsrud et al. 2013; Lind et al. 2018). The southwestern part of the sea is strongly influenced by Atlantic water masses flowing through the area from southwest to northeast. The northern part of the Barents Sea is influenced by cold, Arctic water originating from the Arctic Ocean and from the Kara Sea, and is ice-covered during winter. A coastal current flows along the Norwegian coast into the Barents Sea bringing waters with similar or higher temperature but with lower salinity than the Atlantic water. The “polar front”, a frontal zone with a steep temperature gradient is formed where Atlantic and Arctic water masses meet. The position of the front is mainly determined by topography in the western areas, but is more flexible and less pronounced in the eastern areas. Sea ice covers the areas north of the polar front in winter and may extend to the southern shores in the eastern parts of the Barents Sea.

Fig. 1
figure 1

Map of Barents Sea with topographic and oceanographic features (adapted from Pfirman et al. 1994). Arrows with perforated edges represents sub-surface currents

Full size image

Short description of studied fish stock

Polar cod is ecologically important in the Barents Sea and during the early 1970s, the polar cod was also an important species for commercial fisheries (Ajiad et al. 2011; Hop and Gjøsæter 2013). It is a semipelagic gadid species, often found at near-bottom depths in most of the Barents Sea, excluding the southwestern areas mostly influenced by Atlantic water. While the polar cod is found over large areas in the Barents Sea during the feeding season, spawning is apparently restricted to two separate regions; one in the Pechora Sea in the southeastern Barents Sea (Ponomarenko 1968) and one in the vicinity of Svalbard (Gjøsæter 1973; Boitsov et al. 2013). Recent particle drift studies, (Eriksen et al. 2019; Huserbråten et al. 2019) have corroborated the existence of two separate spawning areas in the region (Fig. 2).

Fig. 2
figure 2

adapted from Huserbråten et al. 2019)

Main spawning areas and drift routes of polar cod eggs and larvae in the Barents Sea, overlaid on the probability of 0-group occurrence from observations (

Full size image

The polar cod is a true arctic species, well adapted to life in ice-covered water at sub-zero temperatures (Rass 1968; Drost et al. 2016; Drost 2017). In the Barents Sea, the early life history of polar cod is particularly associated with ice (Hop and Gjøsæter 2013). Polar cod aggregate in those parts of the Barents Sea where ice forms during late autumn, and spawn at near-bottom depths in areas near the ice edge in late winter-early spring (Ponomarenko 1968). Based on first principles, the buoyant eggs spawned in deeper water will gradually rise towards the surface and aggregate under the ice, possibly in a sub-surface meltwater microlayer. Ponomarenko (1968) stated that since polar cod is the only planktivorous fish found in those areas, and the production of phyto- and zooplankton is particularly rich there, the polar cod lives in a habitat rich in food from its earliest stages of development. In addition, the drift of polar cod arrested under the ice floes is believed to be an environmental adaptation and is comparable to the drift of eggs of other species with currents (Ponomarenko 1964).

The eggs of polar cod are among the largest gadoid eggs (1.53–1.90 mm diameter), which according to Rass (1968), is an adaptation to the cold temperature where the eggs develop. They are also characterized by the egg membrane being very thin, flimsy, and easily damaged. Consequently, they benefit from developing under a protective ice cover where there is no mechanical stress from waves. Rass (1968), also lists late development and scarcity of pigmentation of the embryo to be in accordance with the idea that the ice cover provides protection from direct sunlight.

Material and methods

Description of the data on stock abundance and distribution

Sporadic information about the distribution of polar cod in the Barents Sea based on observations from fishing vessels and from scientific surveys dates to the first half of the previous century (e.g. papers by Andriashev from the 1930s cited by Boitsov et al. 2013). Quantitative information from surveys commenced in the late 1960s to early 1970s, when the acoustic method was developed and was gradually adopted in the Barents Sea. However, the acoustic surveys from the early period were mostly exploratory in nature, and only since 1986 can stock size estimates from annual surveys be compared to recent information (Gjøsæter 2011). Polar cod were normally not the target species for the surveys, which probably did not cover the total distribution area of polar cod, and consequently, the stock size estimates of polar cod should be treated as a relative index of stock abundance (Gjøsæter 1995). The surveys (see description above), in which polar cod were covered, have been conducted in autumn (mostly September) from 1972 to the present, but have undergone major changes over the years. In the early period, the survey was a capelin and polar cod survey, then a pelagic fish survey, then a multispecies survey, and finally, from 2004 onwards, an “ecosystem survey” (Gjøsæter 2011; Eriksen et al. 2018).

Sampling protocol and survey description can be found in Toresen et al. (1998), Michalsen et al. (2013) and Eriksen et al. (2018). In short, the amount of backscatter sampled by the echosounders along the survey track is partitioned among species based on acoustic properties and the species composition in trawl hauls taken along the track. The resulting area backscatter for each species is combined with biological information from targeted trawl hauls (length, weight, age, etc.) to construct a spatial disaggregated stock size estimate distributed on length groups and age groups. The time series from 1986 to 2017 was used here (Fig. 3). For the most part, data over the total distribution area of polar cod were aggregated to form a stock size estimate for the total Barents Sea stock. In addition to the acoustic abundance estimates of polar cod age 1 (1-group) and older, an abundance index of age 0 (0-group) polar cod is constructed from dedicated trawl hauls with a small pelagic trawl towed in the upper 60 m (Eriksen et al. 2009). Prior to 2004, this index was based on a 0-group survey conducted in late August, but since 2004 the 0-group investigation became part of BESS (Eriksen et al. 2018). The timing of the 0-group investigations did not change much.

Fig. 3
figure 3

Abundance indices of polar cod at various life stages, taken from annual survey reports. 0-group indices (eastern and western component separately) and 1-group index (relative numbers) on left vertical axis and acoustic estimate of spawning stock biomass (SSB) on right vertical axis. 0-group and 1-group indices and the corresponding SSB are plotted by year class

Full size image

Description of the modelling activity

The physical habitat experienced by eggs and larvae was derived from a previous particle drift study that described the most likely spawning areas of polar cod in the Barents Sea (Huserbråten et al. 2019). In the prequel study, we applied a Lagrangian particle advection tracking scheme coupled with a 3D dynamical ocean model that logged the physical habitat of eggs/larvae from spawning in spring to 0-group observations in fall. More specifically, particles were released on a regular grid (≈ 40 km equidistance, 537 positions in total) across the entire Barents Sea shelf shallower than 400 m that had been covered by an ice concentration of more than 15% in the period 1990–2017. A new ensemble of 100 particles were released at every point in the grid, every day from 1 January to 30 April, repeated for every year between 1990 and 2017 (yielding a total of 639,030 particles each year). To model the advection of particles in the horizontal plane, we applied the fourth order Runge–Kutta scheme LADIM (Ådlandsvik and Sundby 1994, as for example applied in Myksvoll et al. 2018). As early life stages of polar cod are usually found near the surface (Bouchard et al. 2017), particles were kept at a fixed depth between 0 and 10 m throughout the simulations. Subsequently, an objective search algorithm identified and counted drift trajectories that intersected the 0-group observations of the autumn survey within a three-week period of the surveys and a radius of 20 km. The mechanistic interpretation of the scores from the objective search algorithm (viz. recruitment success to the 0-group) is twofold: first, the backtracking from observations indicates spawning at a given release point for a given year (i.e. supply of recruits from a given spawning location); and second, the frequency of links between spawning location and the observed 0-group abundance may imply a relative survival across the drift phase from the different spawning locations (e.g. see Langangen et al. 2014). For the purpose of the present study, we used the variables temperature and ice concentration above eggs/larvae released from the centre of the most likely spawning area for a given year for the two respective spawning grounds (i.e. the two release locations that had the highest correlation with the 0-group observations), on the day that had the highest probability of reaching the 0-group observations in fall. These temperature/ice profiles thus represented daily averages of 100 eggs/larvae, either spawned east of Svalbard or in the Pechora Sea, over 27 consecutive spawning seasons (from 1990 to 2017) (Fig. 2). To determine the likely time and place of spawning, we identified the day of release that had the highest probability of ending up in the vicinity of 0-group observations in fall. To determine the timing of hatching, we designed an individual-based model (IBM) of egg incubation time as a function of the accumulation of degree days, as shown in laboratory studies by Kent et al. (2016). To test whether timing of spawning, hatching, and ice breakup (i.e. time when ice concentration above eggs/larvae dropped below 50%) differed between the western and eastern spawning areas, we performed two-sample t-tests.

To quantify the relationship between the simulated recruitment success under different environmental conditions and potential predictors, we used a linear regression to model recruitment as a function of Barents Sea ice cover (area of the Barents Sea covered by ice concentration higher than 15%, extracted from a Regional Ocean Modelling System (ROMS) hydrodynamic model as described below), maximum temperature encountered by larvae (10-day mean-filtered over 100 larvae released from a given spawning area), spawning stock biomass (SSB) (estimated biomass of fish above 13 cm was used as an estimate for SSB, Fig. 3), and timing of ice breakup. This regression model was fitted separately for the western (Svalbard) and eastern (Pechora Sea) spawning areas. In the model selection phase, we applied a stepwise model selection scheme with the initial inclusion of all mentioned variables, where only the variables deemed significant were included in the final model. Due to the high degree of collinearity of the ice cover and max temperature covariates, we also did a variation partitioning analysis (Borcard et al. 1992).

The hydrodynamic model used to represent the currents and oceanographic conditions (i.e. temperature, salinity and ice concentration/cover) in the study area was based on the ROMS model, a free-surface, hydrostatic, primitive equation ocean general circulation model (Shchepetkin and McWilliams 2005; Haidvogel et al. 2008). The ROMS model was run with a horizontal resolution of 4 km × 4 km in an orthogonal, curvilinear grid covering parts of the North Atlantic and all the Nordic and Barents seas over the time period 1960–2017 (Lien et al. 2013, as applied in for example Lien et al. 2014). The output from ROMS contained velocity fields, ice concentration, temperature, and salinity in 32 terrain following vertical layers, and a temporal resolution of 24 h.

Results

Seasonal development of the physical environment

Eggs from western and eastern spawning sites (cf. Figure 2) generally experienced similar ice concentrations and temperatures throughout the egg stage, although eggs from eastern spawning sites experienced slightly lower and more variable ice concentration during the first month and during early summer (Fig. 4). Over the 27 spawning seasons modelled, the most likely spawning period started on 7 February and lasted to 9 March with median around 20 February, with no significant difference in timing between the two spawning areas (t = 0.44, df = 50.2, p-value = 0.65). All eggs experienced sub-freezing temperatures (between − 1 and − 2 °C) until start of ice breakup (initiated by the increasing temperatures), which started somewhat earlier in the eastern area yet on average was not significantly different between the two areas (medianwest = 16 July (Julian day 196) vs. medianeast = 22 June (Julian day 172), t = 1.73, df = 53.7, p = 0.08). Despite this slight difference in onset of ice breakup, the average time of hatching was not significantly different between the two areas (t = 1.37, df = 50.7, p = 0.17), and pooled together, the timing of hatching was between 28 April (Julian day 117) and 30 May (Julian day 149) with the median occurring on 10 May (Julian day 129). Finally, in the late larval/juvenile phase (August–September), the temperature in the two areas diverged substantially. Maximum temperatures encountered in the east were significantly higher than in the west (5.7 °C and 2.9 °C, respectively) and were on average between 2.1 and 3.6 °C warmer, respectively (t =  − 7.46, df = 47.7, p-value < 0.0001).

Fig. 4
figure 4

Ice concentration above eggs and larvae that successfully reached the juvenile stage (blue and green bands) and temperature (orange and purple) experienced during drift of eastern and western spawned larvae. Solid lines represent the median and entire polygon the 1st and 3rd quartiles. Spawning and hatching interval represent the minimum, median and maximum times

Full size image

Interannual variability in physical environment linked to 0-group recruitment success

While the two spawning areas exhibited a clear, common seasonal development in their physical environments, there was some interannual variation over the 27-year period, yet with no clear temporal trend (Fig. 5). For example, median ice concentrations above eggs was usually around 90% in most years, but could be as low as ≈40% in extreme years, such as 2005 (Fig. 5a, b). Moreover, larvae/pelagic juveniles from eastern spawning sites generally experienced warmer temperature conditions during the late drift phase compared to western spawning area, with extreme years (e.g. 1995) with median temperatures almost reaching 8 °C (Fig. 5d).

Fig. 5
figure 5

Boxplot of ice concentration above eggs during egg phase of western (a) and eastern (b) spawned eggs, and temperature experienced by larvae/juveniles during late drift phase of western (c) and eastern (d) spawned larvae that successfully reached the juvenile stage

Full size image

The linear regression model for the eastern spawning area showed a high ability to explain larval recruitment success, i.e. a high degree of similarity in environmental exposure (Table 1). Together, ice cover, maximum temperature and SSB explained 74% of the variation in larval recruitment success to the 0-group in the east (Adjusted R2: 0.71, F = 23.2 on 3 and 24 DF, p < 0.0001). During the model selection phase, the variable timing of ice breakup was excluded due to lack of explanatory power. Moreover, the variation partitioning analysis revealed that ice cover and maximum temperature were strongly correlated, although with opposite effects (Fig. 6, Table 1). The SSB also had a significant relationship with recruitment success, indicating a small yet significant stock-recruitment relationship; however, the explanatory power was much lower than the ice cover and maximum temperature. In contrast, the same linear regression model fitted to the western spawning area explained only 0.03% of the variation and was not significant (Adjusted R2: − 0.12, F = 0.03 on 3 and 23 DF, p = 0.99), and none of the variables initially included in the regression analysis could explain a significant portion of the variability in larval recruitment success in the west.

Table 1 Effect of covariates on larval survival (scaled between 0 and 1), from linear regression model fitted to the biophysical models for both spawning areas “Eastern spawning” and “western spawning” (see Fig. 6 for partitioning of variance among the covariates)
Full size table
Fig. 6
figure 6

Venn’s diagram representing the variance partitioning among the covariates ice cover, maximum temperature and SSB in explaining the biophysical model’s ability to predict larval recruitment success in the Pechora Sea (i.e. eastern component)

Full size image

Discussion

During the ongoing warming, climatic conditions in the Barents Sea have changed from moderate in the 1990s to warm temperature conditions in the 2000s and record warm in the 2010s. This warming has led to a strong reduction in the area occupied by Arctic water, which is a preferred habitat for polar cod in the Barents Sea (Ajiad et al. 2011). The redistribution of the occupation area of both juveniles and adults has been linked to this reduction (Hop and Gjøsæter 2013). The distribution of juveniles and adults in the summer-autumn season is well known from joint Russian-Norwegian surveys in August-October, but less is known about where and when spawning takes place (Hop and Gjøsæter 2013; Eriksen et al. 2015). Spawning is known to occur in the southeastern Barents Sea (Ponomarenko 1968; Hop and Gjøsæter 2013), but the 0-group have usually been observed both in the eastern and northwestern area. Our recent particle drift studies have shown that polar cod offspring observed in the northwestern Barents Sea could not come from the southeastern spawning areas (Huserbråten et al. 2019). Furthermore, eggs spawned east of the Svalbard archipelago may drift clockwise around the archipelago and thus could end up in almost all the areas where polar cod juveniles have been observed around the archipelago (Eriksen et al. 2019). This supports earlier findings of prespawners and polar cod juveniles around Svalbard (Ponomarenko 1968; Gjøsæter 1973; Korshunova 2012; Boitsov et al. 2013; Eriksen et al. 2015). Furthermore, Huserbråten et al. (2019) also indicated a clear difference supply of recruits from the two spawning areas, where the Pechora Sea spawning assemblage appeared to be the most important spawning area and the main driver in recruitment variability in the Barents Sea. However, since there was no quantification of the number of adult polar cod spawning in the two respective areas, there was no way to disentangle whether the apparent asymmetry in recruitment success would still hold if scaled by area-specific spawning stock biomass. While an asymmetry in supply from the two spawning assemblages appears plausible, variation in recruitment is likely smaller around the Svalbard archipelago because of the more stable ocean climate east of Svalbard than in the Pechora Sea. In any case, both past and recent observations, along with new simulations of polar cod eggs and larval drift support Ponomarenko’s hypothesis (Ponomarenko 1968) about two separate spawning areas of polar cod in the Barents Sea.

Little is known about the actual spawning period from field observations, but a period lasting for the spring months has previously been inferred from observations of larvae and juveniles and from a few observations of eggs in ice-covered waters. Field observations, summarized by Boitsov et al. 2013, indicate large interannual variations in timing of spawning and hatching in the Pechora Sea. The modelled spawning time in both spawning areas varied by 1 month between years, but variations in the eastern spawning area were more frequent and slightly higher. Hatching time in the model was a function of experienced temperature and incubation time at various temperatures from laboratory studies (Kent et al. 2016). The spawning and hatching periods inferred from the modelling do not necessarily represent the actual spawning and hatching periods, because the model that back-tracked these events are only based on observed juveniles later in the year, combined with the environmental factors they most likely have experienced. This means that only those eggs/larvae that have survived until they are observed in autumn are taken into consideration (i.e. excluding the spawning times that did not lead to successful recruitment to the 0-group). This means that the lack of differences in spawning and hatching periods for the western and eastern spawning areas may be an artefact. Certainly, the results found here are contrary to earlier findings, based on the maturity stage of prespawners, which suggests that polar cod from the Svalbard component mature and spawn later (February–April) than polar cod from the Pechora component (December–March, Boitsov et al. 2013). However, a fully developed mechanistic biophysical model enabling prey- and temperature-dependent growth is needed to increase the quality of the predictions.

Polar cod eggs and larvae are not passive particles. The eggs of polar cod are buoyant (Graham and Hop 1995; Ponomarenko 2000) and will thus quickly rise to the sea surface or the underside of the ice after spawning. The vertical distribution of larvae are poorly documented, but according to Bouchard et al. (2017), 59% of the larvae and juveniles 5–35 mm were found in the top 10 m of the water column from early June to early September in the Canadian Arctic. Field investigations in the 1970s and 1980s also indicated that newly hatched larvae distributed in the upper 10 m in temperatures from 1 to 5 °C during April–May in the Pechora Sea (Boitsov et al. 2013). When the larvae grow larger, they likely distribute over larger parts of the water column. Given the lack of detailed information on larval behaviour, we chose to keep the particles confined to the upper 10 m in our experiment. In preparation for the particle tracking experiment reported on in Eriksen et al. (2019), we tested drift patterns of particles released in 10 m, 20 m and 30 m, and concluded that drift pattern and distance drifted varied little for these three depths. When more detailed information about the behaviours of polar cod larvae becomes available, these behaviours should be taken into consideration in future IBMs of polar cod.

Early fish growth and survival depend on the synchronized development of larvae and their plankton food (Cushing and Horwood 1994). A spring phytoplankton bloom in Arctic water develops when the sea ice starts to melt and forms a stable mixed layer, allowing optimal light and nutrient conditions (Sakshaug and Skjoldal 1989; Wassmann et al. 2010; Oziel et al. 2017) and is followed by a zooplankton bloom with a time lag (Melle and Skjoldal 1998; Loeng and Drinkwater 2007). Oziel et al. (2017) indicated that the spring bloom peaked in May and the summer bloom in September and were more intensive in the Barents Sea waters than Atlantic or Arctic waters. Interannual variability of Chl a in the Barents Sea was found to be about 60%, and variability was linked to large variation in the Air-Ice-Ocean physics (Wassmann et al. 2010; Oziel et al. 2017). In such highly dynamic systems, the match/mismatch between hatching of larvae and availability of their prey is crucial for larval survival (Cushing 1990; Hjort 1914).

The date of ice breakup is probably a very important event in the early life history of polar cod (Fortier et al. 2006; Bouchard and Fortier 2008, 2011; Bouchard et al. 2017). The ice cover is beneficial not only by shielding the eggs from direct sunlight, mechanical stress from breaking waves, and predators (Rass 1968), but also inhibits the onset of the phytoplankton spring bloom and production of food for first-feeding polar cod larvae (Bouchard and Fortier 2008, 2011; Bouchard et al. 2017). When the ice breaks up and melts away, the light penetrates deeper and the meltwater and the rising temperatures of the surface water stabilizes the water column, facilitating the onset of the spring bloom. Synchronization of hatching and ice breakup is probably a prerequisite for successful first feeding, or at least a factor that promotes the survival of early larval stages. The process of ice thawing varies in time and space; however, as explained above, our models estimated dates of ice breakup in the areas where survivors were located and may not reflect the total variability in ice breakup in the Barents Sea as a whole.

We hypothesize that for a species like polar cod that spawns under ice, the timing of hatching of larvae and onset of spring primary and secondary production may be synchronized, since both events probably are triggered by ice breakup through increased light penetration, stabilisation of the photic zone, and increased water temperatures. If this holds true, one may further speculate that the presence of ice will more or less automatically ensure that eggs are hatched at the right time, facilitated by a specific adaption of growth rates that ensures larva are prepared to feed when bloom starts; in contrast to the years when there is little or no ice and egg development has no mediator that dampens their development (i.e. the ice), and they hatch out of synchrony with the bloom. In future studies of polar cod recruitment, these and similar hypotheses should be tested.

While the linear regression model for the eastern spawning area showed a high ability to explain recruitment success to the 0-group, the same model for the western spawning area could not explain a significant portion of the variability in recruitment success. A possible reason for this difference is that for the western area, there is much less variation in sea ice extent and temperatures than in the eastern part of the Barents Sea. The area east of Svalbard is less affected by changes in inflow of Atlantic water and seasonal sea ice extent has not varied as much in this area during the period of observations. Consequently there is less contrast in the environmental data, limiting the ability to detect significant effects in a regression analysis.

The maximum temperatures experienced by the larvae differed significantly between the eastern and western areas; those larvae that survived from the eastern areas developed in significantly warmer water than those from the west. Temperature rose fast in both areas upon ice breakup (cf. Figure 4), but they rose faster and reached a higher maximum in the eastern area. There were also large interannual variations (cf. Figure 5c, d). In most years, the temperature did not rise to a level that exceeds the tolerance level of the larvae and juvenile’s thermal habitat. In 1995, larvae from the eastern area experienced median temperatures of about 8 °C and some larvae as high as 11–12 °C which may be lethal (Koenker et al. 2018). That year was also characterized by the lowest abundance index of polar cod juveniles in the whole survey record. There are reasons to assume that a direct, negative effect of temperature is likely when such high temperatures are encountered. Indirect effects may also be important: In general, when temperatures are within an acceptable range for the larvae, higher temperatures promote growth if sufficient food is available and may also promote the production of food for the larvae and juvenile polar cod, since higher temperatures are often associated with earlier ice breakup and a longer production season. Faster growth promotes survival, because the larvae are susceptible to certain predators for shorter periods, and because faster growth means faster development of sensory organs, ability to capture prey, and escape predators (Hunt et al. 2011). Higher temperatures may, however, be detrimental to larvae that are not able to find suitable food in sufficient concentrations, since increased temperature speeds up the metabolism and increases the demand for food to avoid starvation and death. Higher temperatures can also increase predation rates, since this may imply that polar cod become accessible to predators adapted to higher sea temperatures. Further, increased temperatures may favour more southern lipid-poor zooplankton species over lipid-rich northern species. The idea that a rise in temperatures may have both positive and negative effects on polar cod is corroborated by the findings of Bouchard et al. (2017), that warming temperatures would be good at first, but then would become detrimental when thermal tolerances are exceeded. The maximum temperatures attained in our model in late summer (0.6–5.1 °C in the west and 2.4–9.0 °C in the east) indicate that temperatures encountered are usually within tolerable intervals for larval survival (Koenker et al. 2018). At the same time, the maximum temperature encountered had a large negative effect on recruitment success (cf. Table 1), indicating that polar cod from the eastern spawning areas may in some years encounter temperatures limiting their ability to survive.

In addition to the environmental factors, spawning stock biomass (SSB) was included as a predictor in our models of recruitment success. However, SSB explained little of the variability in recruitment success compared to ice cover and maximum temperature. There are some caveats both regarding data quality and methodology here. First, SSB is poorly determined in the stock assessment of polar cod in the Barents Sea, so the quality of the estimates provided are probably low, both with respect to magnitude and variability. Second, we do not know how spawners are partitioned across the eastern and western spawning areas, so the total SSB were used as a proxy for SSB in both models. The low significance of SSB in the models might result either from a real, low dependence of recruitment magnitude on SSB, or from poor data quality. However, there is likely a strong correlation between SSB and number of eggs spawned. A low dependence of spawning success and resulting year class strength on the biomass of the spawning population are seen for numerous fish stocks, including gadoids characterised by high fecundity. This is not surprising, especially for stocks living in a variable environment, where stochastic environmental effects would often override maternal effects determining survival and growth of larvae.

In summary, the interplay between the physical environment and other ecological factors determining survival is complicated and our current modelling framework was not able to identify which factor plays the most important role: the (1) direct temperature effect on larval physiology; or (2) the overall temperature effect on primary/secondary production and/or composition of plankton community. Also, first-feeding larvae are dependent on a very narrow size spectrum of prey, and high prey abundances do not promote feeding success if prey size exceeds the gape size of first-feeding larvae. That does not mean that temperatures are unimportant; certain temperatures may promote the growth of certain plankton species, and for instance a rise in temperatures may shift spawning areas into more or less suitable areas, or instigate the ice breakup at a more suitable or less suitable point in time. Thus, further studies are needed on the complex interplay between environmental/ecological and physiological effects on polar cod recruitment dynamics.

Profound and sudden changes in the environment like those seen in Arctic areas during recent decades, and those that are forecasted in future climate scenarios, is expected to affect subtle balances in the ecosystems that have evolved over time. Not all changes would be expected to have negative effects, but any change will probably add to the uncertainty of how the ecosystem and its components will behave in the future. Ecosystem components living near the edge of their tolerable environment will be more susceptible to environmental perturbations than others. Polar cod in the Barents Sea live at the southern edge of their natural habitat, and increased dominance of Atlantic water and less ice cover is expected to have negative effects on the stock. The general decrease in stock size in recent decades and the increased variability in year class strength may reflect an overall negative response to warming for this true arctic species in the Barents Sea.