Introduction

The polar cod (Boreogadus saida), a small forage fish, dominates the pelagic fish assemblage in Arctic seas (Fortier et al. 2015). It plays a pivotal role in the transfer of energy from zooplankton to Arctic predators, thus changes in its abundance in response to climate change could alter the services provided to northern communities by the pelagic ecosystem (Welch et al. 1992; Tynan and DeMaster 1997; Darnis et al. 2012). The larvae hatch from January to early July and develop in the epipelagic layer (0–100 m) over spring and summer (Bouchard and Fortier 2008, 2011; Geoffroy et al. 2016). By September and October, age-0 juvenile polar cod initiate an ontogenetic downward migration and recruit progressively to the adult stock in the mesopelagic layer as their size increases (Geoffroy et al. 2016). The abundance and biomass of epipelagic age-0 polar cod in August and September, before this migration, is believed to be a predictor of recruitment to the adult population (e.g., polar cod in Bouchard et al. 2017; other gadids in Laurel et al. 2016).

In the Canadian Arctic, the biomass of juvenile polar cod in August and September increases exponentially with earlier ice breakup (< 50% ice cover) and warmer spring–summer sea surface temperatures (Bouchard et al. 2017). Only the late hatchers survive in years of late ice breakup, resulting in fewer and smaller fish in the fall. In years of early ice breakup, early hatchers survive and enjoy a long growth season, which results in more abundant and heavier fish in the fall. Differences in recruitment to the juvenile stage can be large: juvenile polar cod biomass in September can be 11 times greater for an early May ice breakup compared to a late September ice breakup (Bouchard et al. 2017). Juvenile recruitment was strongly correlated to ice breakup date and sea surface temperatures (SST), two correlated abiotic drivers. However, Bouchard et al. (2017) suspected that the survival of early hatchers in years of early ice breakup is enhanced by biotic factors such as an advanced bloom of ice algae and phytoplankton (e.g., Kahru et al. 2011) and the resulting earlier and more intense production of copepod nauplii and copepodites (Fortier et al. 1995; Ringuette et al. 2002; Daase et al. 2013), the preferred prey of age-0 polar cod (Michaud et al. 1996; Bouchard et al. 2016).

In this study, new data are added (18 region-year combinations in August and 12 in September) to further explore the correlation between ice breakup date and the biomass of juvenile polar cod in late summer reported by Bouchard et al. (2017). We also test the prediction that increased juvenile recruitment in years of early ice breakup is correlated to an earlier phytoplankton bloom and the resulting higher availability of zooplankton. Phytoplankton bloom onset date and acoustically estimated zooplankton standing stock in the fall are used as indices of the production of the pelagic ecosystem over the summer months. Pelagic production and the recruitment of juvenile polar cod are contrasted among three provinces of the Canadian Arctic: Southern Beaufort Sea in the Arctic Ocean Basin proper, the shallow Kitikmeot region in the central Archipelago, and Northwest Baffin Bay including the North Water-Lancaster Sound polynya complex.

Materials and methods

Study area

The Canadian sector of the Arctic Ocean extends from the Beaufort Sea in the West to Baffin Bay in the East (Fig. 1). The sector is divided into several regions by the Canadian Ice Service based on sea-ice characteristics and regime (www.ec.gc.ca/glaces-ice). Typically, these regions are covered by ice for most of the year, with ice breakup occurring from May to September or not at all depending on the year (Bouchard et al. 2017; National Snow and Ice Data Center 2018). From 2006 to 2017, 63 hydroacoustic-trawl surveys of variable duration (1 to 30 days) were completed in ten of these regions in August and September (Fig. 1, Table 1) as part of the ArcticNet annual expedition of the research icebreaker CCGS Amundsen to the Canadian Arctic Ocean (2006–2017) and the Fisheries and Oceans Canada surveys aboard the trawler F/V Frosti in the Beaufort Sea (2012–2014). Overall, valid hydroacoustic estimates of zooplankton and juvenile polar cod biomass were obtained for 40 region-year combinations in August and 23 in September (Table 1).

Fig. 1
figure 1

a Limits of the ten Canadian Ice Service regions analysed in the study grouped by the Southern Beaufort Sea, Kitikmeot and Northwest Baffin Bay oceanographic provinces. b Hydroacoustic and ichthyoplankton surveys from 2006 to 2017. Black lines represent hydroacoustic transects. Colored dots are the locations of ichthyoplankton sampling stations with color indicating the number of net samples at each location over the study period

Table 1 Date range of ichthyoplankton and hydroacoustic sampling periods in each region: Southern Beaufort Sea: MS Mackenzie Shelf, AGM Amundsen Gulf Mouth, AG Amundsen Gulf; Kitikmeot: CM Coronation-Maud, LV Larsen Sound-Victoria Strait, MC M’Clintock Channel, PS Peel Sound; and Northwest Baffin Bay: LS Lancaster Sound, NW North Water, WBB West Baffin Bay

The ten regions surveyed span three main oceanographic provinces (Fig. 1). Southern Beaufort Sea in the west comprises the productive Mackenzie Shelf with intermittent upwelling at its northeast edge (Carmack and Kulikov 1998; Carmack et al. 2004); the mouth of the Amundsen Gulf with the large Cape Bathurst polynya (Stirling 1980; Arrigo and van Dijken 2004; Williams and Carmack 2008); and the deep Amundsen Gulf. The Kitikmeot region in the Southern Canadian Arctic Archipelago (CAA) is characterized by shallow (< 220 m) gulfs, sounds and straits (Coronation-Queen Maud gulfs, Larsen Sound-Victoria Strait), and deeper (< 420 m) sounds and channels (Peel Sound, M’Clintock Channel). Northwest Baffin Bay includes the productive North Water and Lancaster Sound polynya complex, and the western Baffin Bay region (Fig. 1). The west–east general circulation carries surface waters from the Beaufort Sea to Baffin Bay through the shallow CAA (Wang et al. 2012). Biologically, polar cod and the large copepods Calanus hyperboreus and Calanus glacialis dominate pelagic biomass in the deep Beaufort Sea and NW Baffin Bay, while benthic fish and smaller copepods characterize the shallow Kitikmeot (Bouchard et al. 2018, Darnis et al. unpublished data). The North Water and Lancaster Sound polynya complex is considered one of the most biologically productive regions of the Arctic Ocean (e.g., Stirling 1980; Barber et al. 2001; Tremblay et al. 2002).

Remote sensing of ice and Chl a

Ice breakup week (IBW) in a given region-year was defined as the week during which ice concentration fell below 50% (Scott and Marshall 2010) using the Canadian Ice Service data (www.ec.gc.ca/glaces-ice).

For each region-year, the mean chlorophyll a concentration from 1 April to 31 August (Chl a) was estimated using Level 3 daily Aqua MODIS remote-sensing data at a 4 km2 resolution (https://oceancolor.gsfc.nasa.gov/cgi/l3). Daily sea-ice concentration was assigned to each Chl a concentration from the nearest most overlapping pixel of the 25-km resolution Defense Meteorological Satellite Program (DMSP) Special Sensor Microwave Imager (SSM/I)-Special Sensor Microwave Imager/Sounder (SSMIS) (https://nsidc.org, 2006–2016) and of the 12.5-km resolution IFREMER-CERSAT (ftp.ifremer.fr/ifremer/cersat 2017) datasets. Pixels of Chl a concentration with ice concentration > 15% were removed from the analysis to avoid possible contamination of the ocean color signal. For a given region-year, the phytoplankton bloom start date (BSD) was defined as the day when the daily mean Chl a concentration exceeded 20% of the maximum daily mean Chl a concentration from 1 April to 31 August (Platt et al. 2009; Marchese et al. 2017).

Hydroacoustic estimates of polar cod and zooplankton

Ichthyoplankton nets and trawls were deployed (Fig. 1b, Table 1) from the surface to depths varying from 10 to 100 m to ascertain the epipelagic fish assemblage and validate the acoustic signals (details in Bouchard et al. 2017). The fresh standard length (SL) of individual age-0 polar cod was measured on the ship and their weight (W) was calculated based on W = 0.0055 (SL)3.19 (Geoffroy et al. 2016). Data from the two ships for a given region-year were pooled.

Hydroacoustic data were recorded continuously along the track of the ships (Fig. 1b, Table 1) with a Simrad EK60® split-beam echosounder at 38 and 120 kHz (nominal beam angle of 7°). The ping rate varied from ~ 1 to 2 s depending on maximum depth, and pulse duration was set to 1024 µs. Power was 2 kW at 38 kHz and 500 W (2006–2011) or 250 W (2012–2017) at 120 kHz. The co-located transducers were calibrated annually using the standard sphere method (Demer et al. 2015). Conductivity–Temperature–Depth (CTD) profiles from the Amundsen (SBE-911 plus®) and the Frosti (SBE-25® and SBE-19 plusV2®) were used to determine sound speed in water (Mackenzie 1981) and the coefficient of sound absorption (Francois and Garrison 1982) for the acoustic analysis. The echograms were all scrutinized to correct bottom detection by the sounder and to discard noise and signals from other deployed instruments. A time-varied threshold (TVT = 20logR + 2αR − 140, where is the range from the transducer) was also added in the 38 and 120 kHz echograms to offset noise amplification at depth by the time-varied gain (e.g., Geoffroy et al. 2016). A minimum (− 90 dB) and a maximum (− 40 dB) volume backscattering strength (Sv; dB re: 1 m−1) threshold was applied on the data at both frequencies (Benoit et al. 2014; Geoffroy et al. 2016).

The difference in mean volume backscattering strength ΔMVBS (dB re: 1 m−1) between 38 and 120 kHz was used to discriminate polar cod from zooplankton (ΔMVBS120-38 in the range − 10 dB to 5 dB, Benoit et al. 2014; Geoffroy et al. 2016). Monthly mean size (SL and W) of polar cod sampled by nets and the nautical area backscattering coefficient (NASC, m2 nmi−2) in echo-integration cells (0.25 nautical mile long by 3 m deep) at 38 kHz were used to estimate age-0 polar cod integrated biomass (mg m−2) from 13.5 m (effective sampling depth of the transducers) to 100 m. Monthly (August and September) mean integrated age-0 polar cod biomass (B) was calculated for each region-year surveyed (Table 1).

By excluding fish in the top 13.5 m of the water column and at the ice-water interface, our acoustical estimates of age-0 polar cod biomass are underestimating total biomass. The proportion of the age-0 population excluded from the estimates is poorly documented, but the bias is assumed constant across years and regions.

A proxy for zooplankton density in the epipelagic layer (13.5–100 m) was calculated using NASC (m2 nmi−2) in echo-integration cells at 120 kHz. To discriminate zooplankton backscatter from that of fish and macro-zooplankton, only cells with ΔMVBS120-38 > 12 dB (Madureira et al. 1993) were kept in the echo-integration.

Monthly mean zooplankton backscatter (Zoo) was calculated for each region-year.

Copepods dominate the zooplankton of Canadian Arctic seas with the large Calanus hyperboreus and C. glacialis herbivores making up 50 to 90% of the zooplankton biomass and small species such as Pseudocalanus spp., Oithona similis, and Triconia borealis prevailing by numbers (Darnis et al. 2008, 2012). We thus assumed that the acoustic signal attributed to zooplankton is dominated by these species.

All hydroacoustic data analyses were performed with Echoview® 7.1.

Statistics

Spearman’s rank correlations were first tested between abiotic and biotic variables and among biotic variables as exploratory analyses. Each relationship was further investigated using simple linear mixed-effects regression models. In models that included Chl a or zooplankton backscatter as independent variables, these variables were ln-transformed based on a visual examination of scatterplots, as well as on the theory of predator–prey interactions (Holling 1959). Region of sampling was included in the models as a random effect. Relationships between zooplankton backscatter or age-0 polar cod biomass and ice breakup week or phytoplankton bloom start date were also evaluated using second-order mixed-effects models. Linear and second-order mixed-effects models were compared with the Akaike information criteria (AIC). Marginal r2 was calculated for each model (Nakagawa and Schielzeth 2013). Dependent variables in the relationships were ln-transformed prior to statistical analyses to achieve approximate homoscedasticity and normality of model residuals.

Average onset date of the phytoplankton bloom, zooplankton backscatter in August, and age-0 polar cod biomass in August were compared among oceanographic provinces with a one-way ANOVA and a Tukey HSD test.

Statistical analyses were conducted with R® version 3.2.3 (R Core Team 2015).

Results

Ice breakup, microalgal bloom, zooplankton and juvenile polar cod biomass

Across the region-years, the onset date of the phytoplankton bloom was significantly influenced (slope = 0.020, p = 0.002) by the ice breakup week (linear mixed-effects model, Fig. 2a; Online resource 1a). Chl a averaged over the period 1 April to 31 August was significantly correlated to ice breakup week (Spearman’s rank correlation, p < 0.001), although the effect was not statistically significant in the linear mixed-effects model (slope = − 0.018, p = 0.108). Earlier ice breakup tended to result in greater surface microalgal biomass over spring–summer (Fig. 2b; Online resource 1b).

Fig. 2
figure 2

a Bloom start day (BSD) and b average surface chlorophyll a concentration from 1 April to 31 August (Chl a) in relation to ice breakup week (IBW) for all region-years (filled symbols: August; open symbols: September; n = 52) over the period of study. Solid lines are the linear mixed-effects regression models and dashed lines are the 95% confidence intervals. Regression equations are given with the p value of the slope and marginal r2. ρ is the Spearman’s rank correlation

Based on the linear mixed-effects model, zooplankton backscatter in the epipelagic layer (13.5–100 m) in August increased with an earlier ice breakup week (slope =  − 0.170, p < 0.001) and an earlier phytoplankton bloom (slope =  − 0.015, p = 0.014) (Fig. 3a, b; Online resource 2a). An early ice breakup in May increased zooplankton backscatter up to 116-fold relative to a late ice breakup in September (Fig. 3a). A second-order mixed-effects model [ln(Zoo) = 1.11–4.92 IBW–2.89 IBW2] yielded a higher coefficient of determination (marginal r2 = 0.51 vs. 0.36) and a lower AIC value (110.8 vs. 118.7) than did the linear mixed-effects model (Fig. 3a), suggesting that zooplankton backscatter was maximum for an ice breakup in early June. Based on a second-order mixed-effects model as well (marginal r2 = 0.23 vs. 0.17; AIC value = 122.5 vs. 124.4), zooplankton backscatter peaked when the surface bloom started in late May (Fig. 3b). Despite much noise in the relationship, zooplankton backscatter in late summer tended to increase initially with spring–summer Chl a (slope = 0.453, p = 0.363), and then to plateau at Chl a > 1 mg m−3 (Fig. 3c).

Fig. 3
figure 3

Mean epipelagic zooplankton backscatter (Zoo) in August (n = 40) in relation to a ice breakup week (IBW), b bloom start day (BSD), and c mean surface chlorophyll a concentration from 1 April to 31 August (Chl a). Solid black lines are the (a, b) linear or c logarithmic mixed-effects regression models and dashed lines are the 95% confidence intervals. Regression equations are given with the p value of the slope and marginal r2. ρ is the Spearman’s rank correlation. Red curves (a, b) are second-order mixed-effects regression models fitted to the data (see Results)

Biomass of age-0 polar cod (mg m−2) in August in the epipelagic layer increased with an earlier ice breakup (slope =  − 0.254, p < 0.001) (Fig. 4a; Online resource 3a). Polar cod biomass in August was ~ 16 times higher for the earliest ice breakup on week 19 (early May) compared to the latest ice breakup on week 36 (early September). As for zooplankton backscatter, the residuals of the linear mixed-effects regression tended to be positive in June and negative before and after, suggesting that juvenile polar cod biomass in August was maximum when the ice broke up in June. A second-order mixed-effects model also yielded a slightly higher coefficient of determination (marginal r2 = 0.47 vs. 0.44) and a lower AIC value (135.6 vs. 142.6) than did the linear mixed-effects model (Fig. 4a). Age-0 polar cod biomass in August increased with zooplankton backscatter (slope = 0.762, p < 0.001) until it reached a plateau at zooplankton backscatter > 5 m2 nmi−2 (Fig. 4b).

Fig. 4
figure 4

Mean epipelagic age-0 polar cod biomass (B) in August (n = 40) in relation to a ice breakup week (IBW), and b mean epipelagic zooplankton backscatter (Zoo) in August. Solid lines are the a linear or b logarithmic mixed-effects regression models and dashed lines are the 95% confidence intervals. Regression equations are given with the p value of the slope and marginal r2. ρ is the Spearman’s rank correlation. Red curve a is a second-order mixed-effects regression model fitted to the data (see Results). Open symbols represent data presented in Bouchard et al. (2017)

The relationships found in August between abiotic and biotic variables and among biotic variables were also detected in September but at lower statistical significance levels (Online Resources 2 to 5).

Bloom onset date, zooplankton standing stock and juvenile polar cod biomass across oceanographic provinces

The onset date of the phytoplankton bloom averaged over years tended to be more variable in regions of the Southern Beaufort Sea than in the regions of the other two oceanographic provinces (Fig. 5a). Bloom start dates were significantly later (Tukey HSD test, p < 0.05) in the Kitikmeot province of the Canadian Archipelago, and showed no statistical difference between the Southern Beaufort Sea and NW Baffin Bay.

Fig. 5
figure 5

Mean values over the study period of a bloom start day, b epipelagic zooplankton backscatter in August, and c epipelagic age-0 polar cod biomass in August, for different regions of the Canadian Arctic grouped by oceanographic provinces. Number of years analysed over the study period is indicated in bars of panel (a). No values were measured in August in M’Clintock Channel. Oceanographic province with significantly different values based on one-way ANOVA and Tukey HSD test (p < 0.05) is represented by an asterisk

Mean acoustic estimates of zooplankton density in August were similar in the deep Southern Beaufort Sea and deep NW Baffin Bay and significantly lower (Tukey HSD test, p < 0.05) in the shallow Kitikmeot (Fig. 5b).

Juvenile polar cod biomass in August averaged over region-years was highest in Southern Beaufort Sea (mean = 252 mg m−2), intermediate in NW Baffin Bay (mean = 139 mg m−2) and lowest in the Kitikmeot (mean = 19 mg m−2) (Fig. 5c).

Discussion

Early ice breakup and summer pelagic production

In the strongly pulsed primary production cycle of Arctic seas, microalgal biomass first develops in spring at the ice-water interface (Horner 1985; Sakshaug and Slagstad 1991; Søreide et al. 2010). Thanks to early removal of the snow cover over the ice or to leads allowing light to penetrate the surface waters, phytoplankton sometimes blooms before the ice breakup (Haecky et al. 1998; Fortier et al. 2002; Arrigo et al. 2012; Assmy et al. 2017). However, the surface phytoplankton bloom typically starts in the weeks following ice breakup and then deepens towards the nitracline to form a sub-surface chlorophyll maximum (SCM) (Martin et al. 2010). With an increasingly late freeze-up allowing wind mixing, a second surface phytoplankton bloom is now observed in autumn in many Arctic seas (Ardyna et al. 2014).

By allowing light penetration, an earlier ice breakup triggers an earlier phytoplankton bloom (e.g., Ringuette et al. 2002; Kahru et al. 2011; Leu et al. 2011), which results in a longer season of production and, despite the possibility of nutrients becoming limiting in some regions (Tremblay et al. 2012), an overall greater primary production over the spring–summer (Arrigo et al. 2008; Arrigo and van Dijken 2015). Satellite observations detect neither ice microalgae, under-ice phytoplankton, nor deep SCMs (Arrigo and van Dijken 2015). Hence, the remote-sensing measurements of surface Chl a used here certainly underestimated the overall microalgal biomass available to zooplankton grazers over spring–summer (e.g., Assmy et al. 2017). Nevertheless, surface Chl a averaged from April to August, our crude index of ecosystem primary production, was correlated to ice breakup week (Fig. 2). Bloom start date, which dictates the duration of the period of food availability to zooplankton grazers, also depended on ice breakup week as reported previously (e.g., Ringuette et al. 2002; Søreide et al. 2010; Leu et al. 2011).

First-feeding polar cod larvae prey on copepod eggs and the naupliar stages of small copepods such as Pseudocalanus spp. and Oithona similis, whereas juveniles shift their diet to the nauplii and early copepodite stages of the larger Calanus spp. (Fortier et al. 1995; Michaud et al. 1996; Bouchard et al. 2016). An early ice breakup in the North Water relative to Barrow Strait led to the earlier production and higher densities of Pseudocalanus spp. and Calanus spp. early copepodite stages (Ringuette et al. 2002). As well, Leu et al. (2011) observed a temporal mismatch between the algal bloom and the growth of the new copepod generation when ice broke up late in Rijpfjorden (Svalbard) in 2008. These authors suggested that very early ice breakups and algal blooms could also disconnect herbivorous zooplankton from its food, leading to lower population levels, as has been also observed at lower latitudes in the Bering Sea (Hunt et al. 2002).

In the present study, zooplankton backscatter in August was more strongly correlated to ice breakup date and phytoplankton bloom onset date than to Chl a concentration (Fig. 3), indicating that the duration of the season of food availability rather than food abundance likely was the primary driver of the late summer biomass of zooplankton. Over a 4-month range in ice breakup and phytoplankton bloom dates (early May to early September), zooplankton backscatter in August generally increased with an earlier breakup and bloom (Fig. 3). Yet, a close inspection of the relationship reveals that maximum backscatter in August was achieved for ice breakups and blooms occurring in late May and June and that, consistent with the prediction of Leu et al. (2011), the few instances of very early breakup and bloom in May were correlated with somewhat lower zooplankton backscatter in August (Fig. 3a, b). Statistically, the adjustment of a second-order mixed-effects model depicting this maximum in zooplankton recruitment for late May–June breakups and blooms explained a larger fraction of the variance (r2) in zooplankton backscatter in August and resulted in a better model (lower AIC) than a linear mixed-effects model (Fig. 3a, b). Maximum zooplankton backscatter would point to June as the threshold over which climate warming and an ever earlier ice breakup would stop benefiting the present zooplankton assemblage in Arctic seas.

Pelagic production and juvenile polar cod recruitment

Statistical relationships linking fish recruitment to some environmental factor often fail with the addition of new observations, in particular if these come from outside the range of conditions for which the initial relationship was established (Frank 1997; Myers 1998; Leggett and Frank 2008). The addition of 18 new data points for August (from 22 to 40) to the relationship reported by Bouchard et al. (2017) confirmed the link between juvenile polar cod biomass in the fall and the date of ice breakup (Fig. 4a). This suggests that the environmental conditions in Canadian Arctic seas have not yet changed enough to modify the forcing of polar cod recruitment by sea-ice dynamics.

Polar cod is often associated with the ice-water interface and sometimes observed within anfractuosities in the sea ice (Lønne and Gulliksen 1989; Gradinger and Bluhm 2004; Søreide et al. 2006; Melnikov and Chernova 2013; David et al. 2016). This raises the possibility that in region-years of late ice breakup, some age-0 polar cod would ascend to the ice-water interface in August–September, escaping detection in the 13.5 + m layer ensonified by our echosounder. The concentration of juvenile polar cod at the ice-water interface or within the top 13.5 m of the water column could explain in part the estimated low biomasses of age-0 fish in years of late breakup. In the central Arctic Ocean, polar cod sampled at the ice–water interface with the Surface and Under Ice Trawl (SUIT) in August and September were age-1 + fish from 52 to 140 mm SL (David et al. 2016). The 0.3 mm mesh of the inside net of the SUIT should have retained age-0 polar cod if any were distributed at the ice-water interface. Other studies also reported that polar cod associated with sea ice were mostly age-1 and age-2 fish (Lønne and Gulliksen 1989; Melnikov and Chernova 2013). While further studies of sea ice as habitat for age-0 polar cod are warranted (Geoffroy et al. 2016), these observations do not support the hypothesis of a significant migration of age-0 polar cod to the ice–water interface in the fall.

The dependence of juvenile polar cod biomass on ice breakup date can be interpreted in two non-exclusive ways (Bouchard et al. 2017). First, polar cod larvae hatch from as early as January to the first weeks of July (Bouchard and Fortier 2008, 2011). Higher recruitment in late summer could solely result from the fact that an early ice breakup provides early hatching larvae with the minimum temperature and feeding conditions to survive and benefit from a long growth season, leading to abundant and large fish in August. Warmer summer temperatures and more abundant food would play no significant role in polar cod survival and recruitment. Second, in addition to minimum conditions for survival, warmer SST and a general increase in summer pelagic production resulting from an earlier bloom would also contribute to maximize juvenile polar cod biomass in late summer in region-years of early ice breakup. It is difficult to tease apart the respective roles of early breakup and high SST as they are highly correlated (Bouchard et al. 2017). In the laboratory, age-0 polar cod achieved better morphometric and lipid conditions at high temperature (5° vs. 2 °C) or high food ration (2–5 vs. 0.5 prey ml−1) (Koenker et al. 2018). The present field study confirms that the zooplankton food available to young polar cod increases with an earlier ice breakup and phytoplankton bloom onset (Fig. 3a, b); and that the biomass of juvenile polar cod in August is limited at low densities of zooplankton (Fig. 4b). The dependence of polar cod recruitment on zooplankton density followed the expected asymptotic curve predicted by theory when predators become saturated beyond some threshold prey concentration (Holling 1959; Cushing and Horwood 1994). Therefore, the emerging proximal mechanism behind the enhancement of juvenile polar cod recruitment with early ice breakup is the maximization of pre-winter sizes, and possibly lipid storage, by exposing early hatchers to higher SST and saturating feeding conditions.

Bouchard et al. (2017) reported weaker relationships between epipelagic juvenile polar cod biomass and ice breakup date or spring–summer SST for acoustic surveys conducted in October relative to August or September. They attributed the seasonal deterioration of the relationships to the downward migration of the juveniles leaving the epipelagic layer (Geoffroy et al. 2016), and emphasized the importance of conducting surveys during the appropriate temporal window (Bouchard et al. 2017). In the present study, all relationships involving the backscatter of zooplankton were weaker in September than in August (Online resource 2 to 5). By September, the Calanus species making up the bulk of zooplankton biomass have initiated or completed their seasonal vertical migration to depths > 100 m (Dawson 1978; Hirche 1997; Darnis and Fortier 2014). Moreover, an earlier ice breakup and phytoplankton bloom accelerate the development of Arctic copepods and hasten their migration to depth (e.g., Ringuette et al. 2002). While the backscatter of larger macro-zooplankton was excluded in our acoustic analyses, zooplankton backscatter recorded at 120 kHz might have included that of other organisms similar in size to large copepods (e.g., small amphipods). Notwithstanding this potential bias, our results suggest that August is the optimal time window for acoustic surveys aiming at capturing the dependence of juvenile polar cod recruitment on the availability of their epipelagic copepod prey.

Extrapolating the ongoing trend in earlier ice breakup in the different regions studied, Bouchard et al. (2017) projected some admittedly unrealistic increases in the biomass of polar cod by mid-century. They concluded that several factors amplified by the ongoing warming of the Arctic would likely limit an eventual proliferation of polar cod, including the reduction in habitat for the ice-associated individuals and the invasion of Arctic seas by competing or predatory subarctic fish species. In the present study, both the backscatter of zooplankton (Fig. 3a) and biomass of juvenile polar cod (Fig. 4a) were maximum for ice breakups occurring from late May to early July and tended to stagnate or decline for earlier ice breakups in May (note the log scales in both Figs. 3, 4). This parallel response suggests that a mismatch between copepods and their food when the ice breaks earlier than June (Leu et al. 2011) could cascade to the recruitment of juvenile polar cod and contribute to limit the population development of this key species in response to climate warming.

Contrasting oceanographic provinces: depth, seabirds, and the recruitment of polar cod

The Southern Beaufort Sea and the NW Baffin Bay oceanographic provinces are characterized by deep regions and the presence of a large recurrent polynyas (respectively the Cape Bathurst polynya and the North Water-Lancaster Sound polynya complex). The North Water-Lancaster Sound polynya complex in particular is considered an oasis for Arctic predators (Stirling 1980; Brown and Nettleship 1981; Heide-Jørgensen et al. 2013). By contrast, shallow depths (< 100 m) and the resulting absence of warm Atlantic Water limit the abundance of large Calanus copepods and adult polar cod in the Kitikmeot (Bouchard et al. 2018; Darnis et al. unpublished data). Over the period covered by the present study (2006–2017), ice breakups in the regions of the Kitikmeot were among the latest and ranged between early July and early September. Unsurprisingly, the average onset date of the phytoplankton bloom was significantly later (Tukey HSD test, p < 0.05) in this province (Fig. 5a).

Similar to the onset date of the bloom, mean zooplankton backscatter in August presented a clear pattern across the three provinces, with high and similar average values in the deep Southern Beaufort Sea and the deep NW Baffin Bay (6.1–9.7 m2 nmi−2 across regions), and low values (< 2 m2 nmi−2) in the shallow Kitikmeot (Fig. 5b). Late ice breakup and phytoplankton bloom, the low abundance of zooplankton prey, and a scarcity of adult polar cod (Bouchard et al. 2018) likely resulted in low juvenile polar cod biomass in the Kitikmeot (Fig. 5c).

In all regions of the Southern Beaufort Sea and the NW Baffin Bay oceanographic provinces, mean zooplankton backscatter exceeded the 5 m2 nmi−2 threshold (Fig. 5b) under which juvenile polar cod recruitment seems limited (Fig. 4b). Interestingly, despite non-limiting prey density in both provinces, generally higher mean and maximum juvenile polar cod biomasses were observed in the Southern Beaufort Sea than in NW Baffin Bay (Fig. 5c). Mostly because of the availability of nesting cliffs (Gaston et al. 2012), piscivorous seabirds including the thick-billed murre (Uria lomvia), northern fulmar (Fulmarus glacialis) and black-legged kittiwake (Rissa tridactyla), are considerably more abundant in the North Water-Lancaster Sound polynya complex of NW Baffin Bay than in the Southern Beaufort Sea (Wong et al. 2014). The availability of zooplankton prey being non-limiting in the two provinces, we suspect that intense avian predation lowered juvenile polar cod biomass in NW Baffin Bay relative to the Southern Beaufort Sea (top-down instead of bottom-up control). While energy transfer from lower trophic levels to fish and marine mammals through polar cod could be similar in the two provinces, we suggest that intense avian predation on juvenile polar cod in the North Water-Lancaster Sound polynya complex increases overall energy transfer to higher trophic levels. Our results quantitatively confirm the reputation of the North Water as a biological hotspot, further justifying the transformation of the Pikialasorsuaq region into an international protected area under Inuit management (see Eegeesiak et al. 2017).