Climate and juvenile recruitment as drivers of Arctic cod (Boreogadus saida) dynamics in two Canadian Arctic seas

The Arctic is warming four times faster than the global average (Chylek et al., 2022; Rantanen et al., 2022). In the Canadian Arctic (the Beaufort Sea, Canadian Arctic Archipelago, Hudson Bay, and Baffin Bay), this warming has resulted in more ice-free days, thinner sea ice, and warmer summer sea surface temperatures (SSTs) over recent decades (Niemi et al., 2019). Because of relatively low biodiversity in the Canadian Arctic, fewer species control the bulk of energy transfer within the food web compared to other marginal Arctic seas, like the Barents Sea, which could make the marine ecosystems of the Canadian Arctic less resilient to climate and environmental changes (Whitehouse et al., 2014; Murphy et al., 2016; Pedro et al., 2023). Yet, the effects of these environmental changes on marine ecosystems vary between regions. For example, due to more ice-free days, increased stratification, and a change in phytoplankton community, primary production increased in the southern Beaufort Sea, but declined in northern Baffin Bay from the late 1990s to the mid-2000s (Bergeron and Tremblay, 2014; Blais et al., 2017). Further, a comparison of seabird diets collected in the 1970s and 1980s and from 2007 to 2009 suggests that long-term changes in sea-ice cover have resulted in an increase in native boreal forage fish and a decline in Arctic forage species in Hudson Bay, but not in the northern Canadian Arctic (Provencher et al., 2012). Similarly in the low-latitude Canadian Arctic, a decline in sea ice is associated with earlier egg laying since the early 1980s as well as lower chick survival for thick-billed murre (Uria lomvia), but in the high Canadian Arctic there has been no change in reproductive timing (Gaston et al., 2005). Such climate change-induced modifications to species composition, population abundances, population dynamics, or geographic distribution can, in some cases, result in unexpected ecosystem regime shifts (Kortsch et al., 2012; Beaugrand et al., 2015; Fisher et al., 2015). Therefore, a better understanding of the environmental drivers of the dynamics of key marine species is critical for predicting and managing the impacts of climate change on Canadian Arctic marine ecosystems.

Arctic cod (Boreogadus saida) is the most abundant forage fish species in the high Arctic (Benoit et al., 2008; Hop and Gjøsæter, 2013; Geoffroy et al., 2016; Bouchard et al., 2018). As the main prey for groundfish, seabirds, and marine mammals it can funnel more than 70% of the energy between zooplankton and higher trophic levels (Welch et al., 1992). Arctic cod is specifically adapted to the low temperatures and ice conditions of the Arctic. In particular, the early life stages of Arctic cod have a narrow thermal tolerance and rely on sea ice to provide habitat and protection from UV rays, wind and currents, and predators, and therefore the early life stages may be more sensitive to climate change than older Arctic cod (Geoffroy et al., 2023).

Throughout the Arctic, the impacts of climate change on Arctic cod have differed by region. For example in the Barents Sea, the abundance and survival of age-0 fish has declined, possibly due to a decline in sea-ice habitat and increasing surface temperatures, and recruitment is predicted to decline further if conditions continue to warm (Eriksen et al., 2015; Huserbråten et al., 2019; Dupont et al., 2021). Additionally, off the west coast of Greenland, partial recruitment failures have been correlated with increased summer SSTs, though glacial meltwater may provide a thermal refuge for larvae and prevent total recruitment failure as temperatures continue to increase (Bouchard et al., 2020). In contrast, in the Canadian Arctic, warmer SSTs and longer ice-free periods during the summer, coupled with an increase in plankton production, have been associated with higher density and larger size of age-0 Arctic cod, representing a 10-fold difference in age-0 biomass at the end of the first summer compared to years typified by colder conditions (Bouchard et al., 2017; LeBlanc et al., 2020). How increased larval recruitment affects the age-1+ Arctic cod populations in the Canadian Arctic, however, remains unknown.

In other fishes, the survival of early life stages can correlate positively or negatively with the number of fish that recruit to the adult population (Martino and Houde, 2004; Houde, 2008; Laurel et al., 2016; Wilson and Laman, 2021). Therefore, large cohort survival events in Arctic cod may result either in greater recruitment into the age-1+ population or lower recruitment due to an increase in density-dependent mortality caused by greater competition for resources once juvenile fish leave the epipelagic layer to overwinter (Beverton and Holt, 1957; Lorenzen and Enberg, 2002; Martino and Houde, 2004; Dupont et al., 2021). In these contexts, quantifying the sign and magnitude of associations between juvenile and adult population sizes within and among Canadian Arctic cod populations is a critical step towards forecasting future population sizes in response to dynamic environmental drivers.

Environmental forcing associated with global-scale climate indices has been shown to affect the recruitment of early life stages of Arctic cod in the Pacific Arctic (Vestfals et al., 2021). There, shifts in currents and wind speeds, due to changes in climate conditions, can lead to either the local retention or the dispersal of Arctic cod larvae (Vestfals et al., 2021). How these changes in the dispersal patterns of age-0 Arctic cod affect the subsequent dynamics of age-1+ Arctic cod remains unclear. For example, when conditions are warmer over the Pacific Ocean (typical of a negative Arctic Oscillation, a positive Arctic Dipole, and a negative Siberian Alaskan index), a higher proportion of larvae that have been spawned in the North Bering and South Chukchi seas are transported northwest, resulting in a greater contribution of young-of-year Arctic cod to the western Beaufort Sea (Vestfals et al., 2021). This great contribution could lead to a larger age-1+ Arctic cod population in the western Beaufort Sea. Alternatively, these same conditions may also mean greater transport of Arctic cod from the western Beaufort Sea to other regions, like the Canadian Beaufort Sea, then leading to a potential decline in the age-1+ Arctic cod population in the western Beaufort Sea. Given these contrasting potential scenarios and the accelerated rate of warming in the Arctic, providing a better understanding of how environmental variations and climate affect Arctic cod recruitment and population dynamics (age-1+) is imperative.

Here we examine whether sea-ice conditions, summer SST, zooplankton abundance, and climate indices experienced during the early life stages (egg, larva, juvenile) affect the recruitment of Arctic cod into the age-1+ population. We predict that an increase in age-0 Arctic cod and conditions that favor age-0 survival will result in an increase in the biomass of the age-1+ population the following year. We also expect that climatological conditions experienced during the early life stages will impact the age-1+ population the following year. To test the effects of inter-annual variation in biological, environmental, and climatological conditions on age-1+ Arctic cod biomass, we relied on a time series of acoustic-trawl surveys conducted between 2003 and 2019 in the Beaufort Sea and Baffin Bay, and on satellite measurements and climate indices to quantify environmental and climatological variables. In an ecosystem context, variations in Arctic cod biomass are anticipated to propagate through the food web and could result in substantial shifts in marine ecosystem functioning (Harwood et al., 2015; Steiner et al., 2021). Thus, understanding these mechanisms is critical for forecasting Arctic cod dynamics under anticipated climate scenarios.

Acoustic-trawl surveys were conducted aboard the research icebreaker CCGS Amundsen and the F/V Frosti from approximately mid-June through the end of October from 2003 through 2019 (Table 1 and Figure 1). These surveys covered three Canadian Ice Service zone areas of Baffin Bay (North Water, Lancaster Sound, and Northwest Baffin Bay) and the Beaufort Sea (Beaufort Mackenzie Shelf, Amundsen Gulf Mouth, and Amundsen Gulf; Figure 1). The Canadian Ice Service zone areas provide predefined and objective boundaries for survey areas, and make biomass estimates comparable between previous studies that have used the same boundaries (Bouchard et al., 2017; LeBlanc et al., 2020). Survey design was similar to past surveys assessing age-0 Arctic cod (Bouchard et al., 2017; LeBlanc et al., 2020). Acoustic data were recorded continuously by multifrequency Simrad EK60^® or EK80^® echosounders with split-beam hull-mounted transducers at 38, 120, and 200 kHz, all with a 7° nominal beam width. For this study, only data collected using the 38 and 120 kHz transducers were used due to the limited range at 200 kHz. Echosounders were calibrated annually using the standard sphere method (Demer et al., 2015). Narrowband pulses were transmitted and ping rate varied from approximately 1–2 s depending on bottom depth. Pulse duration was set to 1024 µs, and power was set to 2 kW at 38 kHz, 500 W (2006–2011) or 250 W (2012–2019) at 120 kHz. Conductivity-temperature-depth profiles (using the SBE-911 plus^® on the Amundsen and SBE-25 and SBE-19 plusV2^® on the Frosti) were taken during the acoustic surveys to calculate the speed of sound in water and the coefficients of absorption required for acoustic analyses (Simmonds and MacLennan, 2005).

In the Canadian Arctic, Arctic cod generally segregate by size, with age-0 fish in the upper 100 m of the water column and age-1+ fish typically below 200 m (Figure S1; Sekerak, 1982; Benoit et al., 2014; Kjellerup et al., 2015; Geoffroy et al., 2016). This age-1+ population is made up of both immature and mature fish (Geoffroy et al., 2023); for the purposes of this paper this community will be referred to as the age-1+ population. To sample the age-0 fish community and validate acoustic signals, ichthyoplankton and trawl nets were deployed from the Amundsen and Frosti during both day and night during acoustic transects. A double-square net with two square-conical nets (DSN, 1 m² aperture, 500 µm and 750 µm meshes or 1600 µm and 750 µm meshes) and a rectangular midwater trawl (RMT; 8 m² aperture, 1600 µm mesh) were deployed from the Amundsen to sample the ichthyoplankton community. The Frosti did not have a consistent time series of similar ichthyoplankton nets but did deploy a Cosmos-Swan 260 mid-water otter trawl (MWT; headrope = 41.4 m, 12.7 mm mesh codend, Thyborøn Type 2 Standard 5.12 m² steel doors). Therefore, ichthyoplankton that were captured in the MWT when it was towed in the epipelagic (≤100 m) were also used to describe the age-0 community. Similar to previous studies in the region (Geoffroy et al., 2016; Bouchard et al., 2017; LeBlanc et al., 2020), Arctic cod represented all of the age-0 fish with swim bladders caught per year in the epipelagic in both the Beaufort Sea and Baffin Bay (Figures S2a and S3a). Other fish caught in the ichthyoplankton nets lacked a swim bladder, and fish without swim bladders produce weak backscattering signal. Therefore, acoustic backscatter from fish was attributed to Arctic cod.

For biomass estimates of age-0 Arctic cod the average standard length (SL) was calculated monthly by year and by ice-zone area using fish caught with the DSN, RMT, and MWT (Bouchard et al., 2017; LeBlanc et al., 2020). For times (month, ice-zone area, and year combinations) when age-0 fish were not collected, the average SL by ice-zone area and month was calculated using data from all other years when fish were collected during the missing month and in the missing ice-zone area so that there were no missing month-year-ice-zone area combinations.

To examine the age-1+ community composition and validate the source of acoustic signal, trawl nets were deployed both day and night from the Amundsen and Frosti during acoustic transects. An Isaacs-Kidd Midwater Trawl (IKMT) and a beam trawl were deployed from the Amundsen to sample the mesopelagic and demersal fish communities, respectively. The IKMT from 2014–2015 had a 4.5 m × 3.5 m opening, with a 5 mm mesh cod end. A modified IKMT was used in 2016, which had a 3 m × 3 m opening and 12.7 mm mesh cod end. In 2019, the net was replaced with one that had the same mesh configuration but a 4.5 m × 3.5 m opening. The beam trawl on the Amundsen had a 4.3 m footrope and a net with a 9.5 mm codend mesh until 2018, when the net was replaced and had a 6.3 mm codend mesh. On the Frosti, the MWT was used to characterize the mesopelagic fish community. A modified Atlantic Western IIA trawl (AWT; headrope = 22.86 m, footrope = 21.23 m, 12.7 mm codend mesh, Thyborøn Type 2 Standard 5.12 m² steel doors) was deployed to sample the benthic fish communities. Only fish caught by the Frosti that were >4.5 cm, the smallest size collected by trawls on the Amundsen in the Beaufort Sea, were considered for analysis. At least 90% of the fish with swim bladders caught in the mesopelagic and demersal nets per year in the Beaufort Sea and Baffin Bay were Arctic cod (Figures S2b and S3b). Therefore, acoustic backscatter from fish below 100 m was again attributed to Arctic cod.

To calculate biomass, the average size of age-1+ fish by ice-zone area was calculated using fish collected in the beam trawl, IKMT, MWT, and AWT. These nets were only deployed from 2012 onwards. Although not ideal, given a limited time series of data and low sample size in some years, especially in Baffin Bay, the average SL for each ice-zone area was averaged over the period 2012–2019 and applied to the entire acoustic time series, 2003–2019, rather than yearly (Figure 2). We consider the bias from this approach to be negligible because the average yearly SL for each ice-zone area was similar to the combined average for each area (mean difference <1.5 cm between average SL estimates, max difference = 3.9 cm; Figure 2).

Arctic cod biomass and a proxy for zooplankton density were estimated for each year and each ice-zone area in Baffin Bay and the Beaufort Sea. Acoustic data were processed using Echoview^® (Echoview Software Pty Ltd). Methods for analyzing the acoustic data were similar to Bouchard et al. (2017) and LeBlanc et al. (2019). In brief, depending on the type of noise present in the data, denoising algorithms were applied to remove background, transient, impulse, and attenuated noise (De Robertis and Higginbottom, 2007; Ryan et al., 2015) and echograms were visually assessed to manually remove other extraneous noise. The top 18 m of the water column was excluded from analyses to avoid near-field and surface noise. How much of the age-0 biomass is excluded from the analysis is unknown, but excluding the top layer can be assumed to underestimate the biomass of age-0 fish. However, this reduction in biomass should be unbiased and constant (Bouchard et al., 2017; LeBlanc et al., 2020). Backscatter within 1 m of the detected bottom was also excluded from analyses to remove the acoustic dead zone, where distinguishing between biological noise and the bottom echo is difficult (Ona and Mitson, 1996).

Acoustic data were integrated into echo-integration cells 0.25 nmi long by 3 m deep. A maximum (−40 dB) volume backscattering strength threshold was applied to the echo-integration cells. This threshold is greater than the densest schools of Arctic cod and most likely is an artifact of non-biological noise (Geoffroy et al., 2016). A minimum threshold (−80 dB) was applied to the 38 kHz data below 150 m to exclude backscatter from zooplankton. Therefore, the majority of the 38 kHz signal below 150 m can be assumed to have originated from fish (Geoffroy et al., 2016). Above 150 m (the effective sampling depth of the 120 kHz transducer due to low signal-to-noise values during some years), cells with a difference in mean volume backscattering strength (ΔMVBS; dB re 1 m⁻¹) at 120 kHz and 38 kHz that was greater than 5 dB were assigned to zooplankton, and values between −10 dB and 5 dB were assigned to fishes with swimbladders (Benoit et al., 2014; Geoffroy et al., 2016; Bouchard et al., 2017; LeBlanc et al., 2020).

The weight-based target strength function (TS_W; dB re 1 m² g⁻¹) was used to estimate age-0 and age-1+ Arctic cod biomass (Simmonds and MacLennan, 2005; Benoit et al., 2008):

where TS_N is the estimated target strength (dB re 1 m²) and $W_$ is the average weight (g).

For each age group (age-0 and age-1+), the TS_N was estimated based on the average Arctic cod SL. TS_N was calculated using the following equation from Geoffroy et al. (2016):

where TS_N is the target strength per individual (dB re 1 m²) and $SL_$ is the mean SL (cm) of each age group.

The estimated average weight (⁠$W_$⁠) was calculated using the following equation used in Benoit et al. (2008):

where SL_j is the average standard length of length class j contributing the fraction n_j of the total length-frequency distribution of the population (Figures S4–S7) and ΔSL (1 cm) is the difference between length classes. Values for a_f (0.0055) and b_f (3.19) are from the Arctic cod length-weight regression described in Geoffroy et al. (2016).

Biomass (B; g m⁻²) was then estimated for each echo-integration cell using the nautical area backscattering coefficient (NASC, m² nmi⁻²) at 38 kHz and TS_w using the following equation:

Age-1+ fish biomass was then integrated (g m⁻²) across the survey period over depths between 100 m and 550 m, the deepest effective sampling depth of the 38 kHz transducer. Although age-1+ fish biomass is underestimated by excluding depths below 550 m, peak biomass typically occurs between 300 and 500 m across the Canadian Arctic (Majewski et al., 2016; Majewski et al., 2017; LeBlanc et al., 2019). Integrated biomass was also calculated for age-0 fish down to 100 m from July through September. Average biomass (age-0 and age-1+) was then calculated for each year and ice-zone area to compare yearly estimates for each region. The integrated NASC (m² nmi⁻²) at 120 kHz, obtained from acoustic data assigned to zooplankton based on the ΔMVBS method previously described, from July through September in the depth range of 18–150 m was calculated as a proxy for zooplankton density (LeBlanc et al., 2020).

Similar to Bouchard et al. (2017) and LeBlanc et al. (2020), environmental factors included ice-breakup week and average summer SST for each year and ice-zone area. Using data from the Canadian Ice Service database (Ice Graph, 2021), ice-breakup week was defined in each area as the first week when total ice concentration reached less than 50% (Bouchard et al., 2017; LeBlanc et al., 2020). Average summer SST was calculated from May 1 to July 31 (Bouchard et al., 2017; LeBlanc et al., 2020) using level-4 satellite-derived estimates of daily SSTs at 0.05° × 0.05° resolution for each of the ice-zone areas from the Copernicus website (Copernicus Climate Change Service, 2022).

The Arctic Dipole (AD), Arctic Oscillation (AO), and North Atlantic Oscillation (NAO) indices were selected to determine if climate indices during the early life stages influence the recruitment of age-1+ Arctic cod. The AO and NAO are two of the major drivers of atmospheric circulation and pressure in the Arctic (Hanna et al., 2021). The NAO is the fluctuation of pressure of the North Atlantic and is the difference in pressure at sea level (SLP) between the Icelandic Low and Azores High, while the AO is the fluctuation in pressure over the Arctic and is the difference in pressure between the Arctic and mid-latitudes. Often the effects of the NAO and AO on the Arctic are the same (Dickson et al., 2000), and fluctuations in the NAO and AO are linked to changing temperatures, precipitation patterns, and wind direction and strength. The AD has recently been recognized as a third important climate index in the Arctic, influencing meridional wind patterns (north–south) rather than zonal wind patterns (east–west) like the NAO and AO (Overland and Wang, 2010; Overland et al., 2012). The AD index was calculated as the first empirical orthogonal function (EOF) pattern of 70–90°N regressed to the SLP anomaly (Wu et al., 2006); data were downloaded from the NOAA Bering Climate website (Bering Climate, 2021). The AO index was calculated as the EOF pattern of SLP from 20–90°N regressed to the SLP anomaly. The NAO index was calculated as the first EOF pattern of SLP from 0–90°N regressed to the SLP anomaly. Both NAO and AO data were downloaded using the rsoi R package (Albers and Campitelli, 2019) which facilitates the retrieval of data from the NOAA Climate Prediction website (Climate Prediction Center, 2022). Average climate indices were calculated for each season: Winter (December, January, February), spring (March, April, May), summer (June, July, August), and fall (September, October, November). Average climate indices were also calculated for each year.

Generalized additive models (GAM) were used to analyze the lagged effects (those that would have influenced the early life stages) of changing environmental, biological, and climatological conditions on age-1+ Arctic cod biomass. GAMs are flexible extensions of generalized linear models that model the relationship between the response variable and predictor variables as nonparametric smooths (splines; Wood, 2004; 2006). GAMs were fit for both regions (Beaufort Sea and Baffin Bay) and for each type of variable: environmental (SST, ice-break-up week), biological (age-0 biomass and zooplankton NASC), and climatological (AD, AO, and NAO) (Table S1). The effects of environmental variables were modelled separately due to their strong correlations with each other in each region (r_s > −0.7). The effects of the climate indices were also modelled separately to avoid issues of co-linearity between years and having too many covariates with limited sample size. Variables were tested with a lag of 1 year and 2 years to represent the conditions experienced during the early life stages of Arctic cod. Variables were only lagged 1 year or 2 years for three reasons: (1) Arctic cod is short-lived (max age = 7; Hop et al., 1997), (2) environmental factors have been shown to have the strongest effects on ages one and two (Dupont et al., 2021), and (3) sample size was limited when biological variables were lagged by more than 2 years. Average age-1+ biomass data were ln-transformed and models were fit with a Gaussian distribution with an identity link function. In each model, year was also included as a random effect to account for yearly variation not accounted for by other variables (e.g., number of survey days or distance surveyed).

The environmental and biological effects were modelled with thin plate regression splines, and climatological effects were modelled with cubic regression splines. To reduce potential overfitting and to increase interpretability the number of knots was limited to five for each spline. Because model likelihoods with different fixed effects cannot be compared when fit using restricted maximum likelihood (REML) methods (Zuur et al., 2009), models were initially fit with generalized cross validation (GCV), which is a measure of the model prediction error and a consistent means of model selection (Hughes et al., 2015). Climatological and environmental models were compared and model selection was based on Akaike information criterion (AIC) scores, GCV scores, and the adjusted r² value. Models with the lowest AIC and GCV scores, and highest adjusted r² values, were chosen. Biological models could not be compared to one another because not all ice zones were sampled each year, and so the model fit with biological data lagged by 1 year had a different sample size than the model lagged by 2 years. For example, in Baffin Bay biological models with variables lagged by 1 year had a sample size of 25, while models with variables lagged by 2 years had a sample size of 23, and environmental and climatological models both had a sample size of 36.

Best-fitting models were re-estimated using REML methods, which is more robust to under-smoothing (Wood, 2006). Variable selection was based on the shrinkage method using a double penalty approach, which can shrink all functions towards zero, essentially removing unimportant variables from the model and leaving a more parsimonious model (Marra and Wood, 2011). Variables that were shrunk towards zero were removed from the final model. Model residuals were checked for first-order autocorrelation using the autocorrelation function. All analyses were conducted in R 3.6.3 (R Core Team, 2019) with the mgcv package (Wood, 2015). Figures were created with the ggplot2 (Wickham, 2016), ggmap (Kahle and Wickham, 2013), ggpubr (Kassambara, 2020), and cowplot (Wilke, 2019) packages.

Mean age-1+ Arctic cod biomass varied through time and by region (Figure 3). Generally, mean age-1+ biomass (averaged over area-years) was higher in Baffin Bay (mean = 3.6 g m⁻²) than in the Beaufort Sea (mean = 1.13 g m⁻²). However, the variance around average biomass estimates for each area was greater in Baffin Bay than the Beaufort Sea (Figure 3), which could be due to the lower transect coverage in Baffin Bay (Table 1). In the Beaufort Sea the average age-1+ biomass was lowest in the Mackenzie Shelf area (mean = 0.46 g m⁻²), and in Baffin Bay it was lowest in the West Baffin Bay area (mean = 2.66 g m⁻²). In the Beaufort Sea mean age-1+ biomass (averaged over area-years) was greatest in 2011 (mean = 3.3 g m⁻²) and 2017 (mean = 3.2 g m⁻²). In Baffin Bay mean age-1+ biomass (averaged over area-years) was greatest in 2013 (mean = 9.0 g m⁻²; Figure 3). Inter-annual variation in average biomass estimates was high, with a median percent change between years of 97% in the Beaufort Sea and 94% in Baffin Bay. Although age-1+ Arctic cod biomass was greater in Baffin Bay, age-0 biomass was approximately 2.3 times greater (averaged over area-years) in the Beaufort Sea than in Baffin Bay (Figure 4). In Baffin Bay, age-0 biomass was lowest (averaged over area-years) in the West Baffin Bay area (mean = 0.10 g m⁻²), and in the Beaufort Sea, it was lowest (averaged over area-years) in the Mackenzie Shelf area (mean = 0.31 g m⁻²).

In the Beaufort Sea the best-fit GAMs explained 57%–87% of the deviance (Tables 2, S1, and S2). The best-fit biological model included biological variables lagged by 1 year (age-0 biomass and zooplankton NASC), explained 87% of the deviance, and revealed that age-0 biomass lagged by 1 year had a significant positive effect (F = 105, p = 0.003) on age-1+ biomass, while the backscatter from zooplankton did not have a large effect (F = 2.8, p = 0.189; Tables 2 and S1 and Figure 5a). In comparison, the model with biological variables lagged by 2 years only explained 55% of the deviance (Table S2), and there was a weak negative but non-significant relationship (F = 2.79, p = 0.071) between age-0 biomass lagged by 2 years and age-1+ biomass (Figure S8). All climate models had similar model diagnostics and fits (Table S2), possibly due to the underlying relationship between the NAO and AO. However, the model with seasonal AO indices lagged by 1 year had the best fit and explained 70% of the deviance. Seasonal AO indices (spring, summer, fall) lagged by 1 year had significant effects on age-1+ biomass (Table 2, Figure 5b–d), with the AO index in spring (F = 9.94, p < 0.0001) and in fall (F = 9.88, p = 0.0003) having the greatest effects. Age-1+ biomass was lowest when spring AO was between 0 and 1 (Figure 5b). Age-1+ biomass decreased as summer AO increased until AO values were greater than zero, when biomass started to increase again (Figure 5c). Age-1+ biomass decreased as fall AO indices increased (Figure 5d). The best fit GAM that included environmental variables (ice-breakup week and SST) included SST lagged by 1 year; however, the lagged SST had a weak and non-significant effect (F = 0.098, p = 0.432) on the age-1+ biomass (Table S2).

For Baffin Bay, the best-fitting GAM models explained 48%–71% of the deviance (Tables 2, S1, and S2). The GAM model with biological variables lagged by 1 year only explained 48% of the deviance and suggested that the effect of age-0 biomass on age-1+ Arctic cod biomass was weak (F = 0.152, p = 0.248; Table S2). The best-fit GAM that included biological variables lagged by 2 years (deviance explained = 65%) and revealed that the zooplankton backscatter had a non-significant negative (F = 2.22, p = 0.067) effect on the age-1+ biomass (Table S2, Figure S9); however, sample size was small (n = 23). Like in the Beaufort Sea, model diagnostics were similar between climatological models (Table S1), but the model with seasonal NAO values lagged by 1 year gave a better fit, explaining 57% of the deviance in age-1+ Arctic cod biomass. Similar to the Beaufort Sea, the spring (F = 5.96, p = 0.005), summer (F = 4.53, p = 0.042), and fall (F = 5.49, p = 0.013) climate indices had significant effects on the age-1+ biomass, with spring and fall having the greatest effects (Table 2 and Figure 6a–c). Estimated smooths suggest that age-1+ biomass was highest when the spring NAO index was either positive or negative, but lowest at zero (Figure 6a). However, age 1+ biomass decreased as summer and fall NAO became more positive (Figure 6b and c). The best fit environmental GAM (explained deviance = 0.71) suggested that there was a significant effect of ice-breakup week (F = 1.24, p = 0.031) on age-1+ biomass, with the highest age-1+ biomass occurring when ice breakup was early (21 weeks) or late (32 weeks), and the lowest biomass when the ice-breakup week was in between (26 weeks; Figure 6d).

The Canadian Arctic food web depends largely on the productivity of Arctic cod to funnel energy from lower levels. This study demonstrates that Arctic cod is tightly linked to environmental conditions and climate, supporting the idea that the community structure and variation in energy flow in the Canadian Arctic could be driven by bottom-up processes (Carmack and Wassmann, 2006; Conservation of Arctic Flora and Fauna, 2013). There have been few studies that attempt to examine age-1+ Arctic cod biomass in the Canadian Arctic over a period longer than a decade, but acoustic surveys conducted aboard the CCGS Amundsen and the F/V Frosti provided a unique and robust dataset, which allowed for an in-depth assessment of trends in Arctic cod biomass over relatively large spatial and temporal scales.

The seasonal NAO and AO during the first year of development had significant effects on the biomass of age-1+ Arctic cod in Baffin Bay and Beaufort Sea, respectively. In addition, the random effect of year was reduced to zero in the models but was retained in all the environmental and biological models, suggesting that the random effect of year in those models may account for some of the variance due to seasonal climate oscillations. Although the effects of the NAO and AO are sometimes indistinguishable from one another (Dickson et al., 2000), the NAO may have a lesser effect on Arctic cod biomass in the Beaufort Sea because the center of the NAO moves and is less influential over the Pacific Arctic outside of winter (Folland et al., 2009; Bednorz et al., 2019). The abundance and biomass of other forage fish species, such as capelin (Mallotus villosus), Pacific herring (Clupea pallasii), and Antarctic silverfish (Pleuragramma antarcticum; Hjermann et al., 2004; Puerta et al., 2019; Murphy et al., 2021; Corso et al., 2022), as well as other gadid species, Pacific cod (Gadus macrocephalus), Atlantic cod (Gadus morhua), and walleye pollock (Gadus chalcogrammus; Puerta et al., 2019; Zimmermann et al., 2019) in polar and sub-polar seas have also been shown to vary with decadal climate indices.

The effects of the NAO are typically more pronounced in winter (Hurrell and Van Loon, 1997), but here the spring and fall indices had the greatest influence on the biomass of age-1+ Arctic cod for both the NAO and AO (Figures 5 and 6 and Table 2). In the Canadian Arctic, peak hatching time for Arctic cod is during the late spring and early summer (Bouchard and Fortier, 2011), after which they remain in the epipelagic layer for several months where they are vulnerable to surface conditions (Benoit et al., 2014; Geoffroy et al., 2016). Positive NAO and AO conditions are often associated with a colder Arctic. In the Beaufort Sea during years when a positive NAO or AO persist past the spring, more ice is pushed into the region (Drobot and Maslanik, 2003; Armitage et al., 2018). These conditions may inhibit the growth and survival of age-0 Arctic cod, ultimately reducing the number that recruit into the age-1+ population. Some late-hatching fish also recruit to sea ice (Geoffroy et al., 2016), and these fish may be pushed further into the region along with the sea ice and into less suitable habitats like the Canadian Arctic Archipelago.

In the Bering and Chukchi seas, variation in the distribution and retention of age-0 Arctic cod has been correlated to regional advection and large-scale climate forcing, most likely due to wind-driven currents and the variation in trajectory and velocity of these currents (Levine et al., 2021; Vestfals et al., 2021). In the Beaufort Sea, during years of strongly positive and negative AO, the difference in sea level pressure generates cyclonic or anticyclonic circulation anomalies, respectively, with winds and along-shelf geostrophic circulation anomalies pushing surface waters and sea ice east (positive AO) or west (negative AO; Proshutinsky et al., 2015; Armitage et al., 2018; Wang and Danilov, 2022), with current speed anomalies increasing 0.5 cm s⁻¹ per AO index (Armitage et al., 2018). Similarly, during years of strongly positive or negative NAO in Baffin Bay, there is stronger or weaker cyclonic circulation, respectively, with more (positive NAO) or less (negative NAO) Arctic waters flowing south into Davis strait (Münchow et al., 2015). We hypothesize that, during years with positive winter or spring NAO and AO, more age-0 fish are carried with the sea ice and currents from the western Beaufort Sea into the eastern Beaufort Sea, like the Amundsen Gulf, and more age-0 fish are carried into northern and western Baffin Bay from eastern Baffin Bay (Figure 7). Thick sea ice may also weaken the effects of wind and current anomalies during winter and spring (Armitage et al., 2018), but if the positive phase continues into summer and fall, when there is also less ice to modulate the effects, the distribution of age-0 fish could shift dramatically, carrying more age-0 fish towards less suitable habitats. For example, with continued transport of fish east in the Beaufort Sea during positive AO years, more age-0 fish may end up in the Canadian Arctic Archipelago where the survival rate is low due to shallow water and thick sea ice (Figure 7; Bouchard et al., 2018). In Baffin Bay, a positive NAO could result in more fish transported out of Baffin Bay into less suitable habitats further south, such as the Labrador Sea (Figure 7).

Otolith microchemistry studies have shown that the early life stages of Arctic cod can drift for long distances (Bouchard et al., 2015), but genetic studies suggest that oceanographic conditions separate western and eastern Arctic cod populations (Nelson et al., 2020). However, fish in eastern Canadian populations, from Resolute Bay to the Gulf of St. Lawrence, were genetically similar (Nelson et al., 2020), which may support the hypothesis of southward transport. In addition, a positive NAO has been correlated positively with a stronger Labrador current (Han et al., 2014), as well as with the number of icebergs off Newfoundland (Drinkwater, 1996; Münchow et al., 2015; Cyr and Galbraith, 2021), and during colder years (indicative of a positive NAO index) there are more Arctic cod caught off the coast of Newfoundland, the most southern limit of their distribution (Lilly et al., 1994; Brown et al., 1996; Rose, 2003). These findings support the idea that during years of a more positive NAO, Arctic cod may be more likely to be advected south from Baffin Bay.

Similar to gadids in other regions (Laurel et al., 2016; Wilson and Laman, 2021), the relationship between the lagged biomass of age-0 Arctic cod and the age-1+ biomass in the Beaufort Sea was positive. This relationship was not evident in Baffin Bay, possibly because age-0 biomass is lower in this region (Figure 4), explaining the lack of relationship between age-0 biomass and age-1+ biomass measured in this study. In addition, the average size of fish in Baffin Bay was longer than in the Beaufort Sea, which could suggest that there are more older fish in Baffin Bay and that the relationship between older fish (age-2+) and age-0 fish may be weaker than that between age-0 fish and age-1 fish. Models for Baffin Bay explained less variance overall (48%–71%), compared to Beaufort Sea models (57%–87%), which could also be due to fish in Baffin Bay potentially being older than in the Beaufort Sea; thus any relationship influencing the early life stages would explain less variance in the Baffin Bay populations. The lower explanatory power of models in Baffin Bay could also be due to higher variance in age-1+ biomass estimates (Figure 3). This higher variance could be due to more spatial variation in age-1+ Arctic cod distributions in Baffin Bay or because there were fewer transects over a larger area in Baffin Bay (Table 1) making trends harder to characterize. It could also be due to Baffin Bay being more complex than the Beaufort Sea, with more inflow of Atlantic water and boreal species from the Atlantic Ocean which could create more chances for competition with other species or higher risk of predation.

The Beaufort Sea hosts greater biomass of age-0 fish, possibly because it provides better nursery habitats for age-0 Arctic cod relative to more eastern areas. Similar to results in LeBlanc et al. (2020), which used a shortened time series of the data included here, the biomass of age-0 fish was approximately 2.3 times greater in the Beaufort Sea (Figure 4) than in Baffin Bay. In the Canadian Arctic, Arctic cod hatch under the ice in the epipelagic layer from January through July, and then remain there until late fall, when they join their age-1+ congeners at depth to overwinter (Bouchard and Fortier, 2011; Geoffroy et al., 2016). However toward the end of their first summer, age-0 fish are highly susceptible to predation (LeBlanc et al., 2019; LeBlanc et al., 2020), which might be much higher in Baffin Bay compared to Beaufort Sea due to higher abundances of top predators such as thick-billed murre and northern fulmar (Fulmarus glacialis) (Wong et al., 2014). In addition, there is a greater influx of freshwater in the Beaufort Sea (Bouchard and Fortier, 2011; Bouchard et al., 2015), with the Mackenzie River feeding into the Mackenzie Shelf, making it the most estuarine shelf in the Canadian Arctic (Macdonald et al., 1987). This increased freshwater supply leads to warmer winter temperatures needed for faster egg development, earlier under-ice hatching, and greater larval motility, all potentially leading to the increased survival of age-0 fish during the open water period (Bouchard and Fortier, 2011). These early hatching fish in the Beaufort Sea further benefit from the warmer summer SSTs (mean = 3.14°C), compared to those in Baffin Bay (mean = 1.45°C; Bouchard et al., 2017; LeBlanc et al., 2020), which may be why age-0 fish are larger in the Beaufort Sea (Figures S4 and S5).

To assess biomass, a demarcation of 100 m was designated as the cutoff of age-0 and age-1+ Arctic cod; however, there are some uncertainties in the biomass estimates for the two groups. For instance, while the majority of age-0 fish remain in the epipelagic layer over the summer and fall and the age-1+ fish remain below 200 m, a small percentage of age-0 fish will descend below 100 m earlier in the year (Sekerak, 1982; Benoit et al., 2014; Kjellerup et al., 2015; Geoffroy et al., 2016), and these fish may get included in the age-1+ biomass estimates. However, the number of age-0 fish is relatively low compared to the number of age-1+ fish (Sekerak, 1982; Benoit et al., 2014; Geoffroy et al., 2016). For example, Sekerak (1982) sampled Lancaster Sound and western Baffin Bay and found that the densest age-0 aggregations (10.4–19.8 fish per 100 m⁻³) occurred at depths of 10–20 m while age-0 fish below 100 m were much less dense (0.38–1.22 fish per 100 m⁻³). Therefore, due to their small size and small numbers, their contribution to the age-1+ biomass estimates should be minimal. Furthermore, due to limited length-weight data from the net surveys, we used a consistent length-weight relationship from Geoffroy et al. (2016). These length-weight relationships could change over time due to prey quality or availability, or due to environmental or climatological conditions. Using a consistent length-weight relationship assumes that the condition of the fish remained the same through time, which could introduce additional uncertainties around the biomass estimates. Other sources of uncertainty in biomass estimates may also be due to the combination of EK60 data from 2003–2016 and EK80 data from 2017–2019 collected on the Frosti in the Beaufort Sea. Due to over-amplification of low-power signal by the EK60, the biomass estimates from the EK80 can be 3%–12% lower than those produced by the EK60 depending on the depth, frequency, and scattering strength of the target (De Robertis et al., 2019). This difference may make comparing biomass estimates between years that used the EK60 (2003–2016) and EK80 (2017–2019) more difficult. However, this issue only affects 3 years in the Beaufort Sea; with the variation around the mean biomass estimates being high, it should incorporate the potential 3%–12% difference.

In contrast to age-0 biomass, the biomass of age-1+ Arctic cod was approximately three times greater in Baffin Bay than in the Beaufort Sea (Figure 3). This difference may be because the relationship between age-0 biomass lagged 2 years and age-1+ biomass in the Beaufort Sea was weakly negative, although not significant (Table S2 and Figure S8), possibly due to higher age-0 biomass causing density-dependent competition resulting in lower recruitment. The Beaufort Sea is less productive (total areal primary productivity = 8 Tg y⁻¹) than Baffin Bay (62 Tg y⁻¹; Sakshaug, 2004). Additionally, regions like the North Water and Lancaster Sound in Baffin Bay are known for being highly productive (Stirling, 1980; Heide-Jørgensen et al., 2013). The Eastern Canadian Arctic may thus provide better resources throughout the year for Arctic cod, resulting in higher age-1+ Arctic cod abundance. However, there was a weak negative correlation between zooplankton lagged 2 years and age-1+ biomass in Baffin Bay, which may suggest that some sort of density-dependent mortality is happening similar to that seen in the Beaufort Sea with the relationship of age-0 fish lagged by 2 years. An increase in zooplankton abundance was correlated previously with an increase in size and biomass of age-0 Arctic cod at the end of their first summer (LeBlanc et al., 2020). However, this relationship did not translate to more age-1+ biomass in this study, which could be due to the fact that we averaged zooplankton backscatter over a longer time period (July through September) than the previous study (LeBlanc et al., 2020), in which monthly zooplankton and larval indices were analyzed for a finer temporal scale resolution.

In Baffin Bay, age-1+ Arctic cod biomass was significantly affected by the timing of ice breakup during early life stages. Previous research, further supported by the inclusion of more data in this study, has shown that earlier ice breakup results in greater age-0 biomass (Figure S10; Bouchard et al., 2017; LeBlanc et al., 2020). Interestingly, however, our study suggests that either early ice breakup or late ice breakup results in higher age-1+ biomass while median breakup timing does not confer any benefits to the age-1+ population. The years with early ice breakup most likely benefit earlier hatching fish and provide a longer growing season for Arctic cod and more zooplankton prey leading to increased survival and recruitment into the age-1+ population (Bouchard et al., 2017; LeBlanc et al., 2020). However, during years of late ice breakup only late-hatching fish survive (LeBlanc et al., 2020). Although this survival effect may lead to smaller and fewer fish at the end of the first summer, it seems to lead to increased age-1+ Arctic cod biomass. These late-hatching fish may be more likely to recruit to sea ice at the end of their first summer, rather than descending to depth to overwinter with their age-1+ congeners, which may extend their growing season and lead to decreased competition and predation through the winter and ultimately increased survival into the age-1+ population (Craig et al., 1982; Geoffroy et al., 2016). By joining their congeners at depth the following year when Arctic cod are larger, their larger size may make them able to better compete for food and avoid predation and contribute more biomass than smaller cod that descended to depth in their first season and were eaten or starved.

In the Beaufort Sea, where there is a strong positive relationship between age-0 biomass and the age-1+ population of the following year (Figure 5a), the warmer conditions that favor age-0 fish may continue to benefit the age-1+ population, at least in the short term. However, in Baffin Bay, where there was not a strong relationship between age-0 and age-1+ biomass, continued warming may not confer the same benefits to the age-1+ population in the region. Ice-breakup week is predicted to happen earlier through time (Stroeve et al., 2012; Bouchard et al., 2017) and could occur during week 26 through 28 in northwest Baffin Bay as predicted by Bouchard et al. (2017). This earlier breakup could further contribute to a decline in age-1+ biomass in that area. However, due to the complex parabola-like relationship between ice-breakup week and age-1+ biomass, with ice-breakup week occurring in the middle having the lowest biomass, some areas in Baffin Bay may eventually see an increase in age-1+ biomass if ice-breakup week occurs as early as week 22 or 23 as predicted for the North Water (Bouchard et al., 2017). However, there was a strong negative effect of a positive AO and NAO (Figures 5b–d and 6a–c) in the Beaufort Sea and Baffin Bay, respectively, and both climate indices are forecast to become more positive (Gillett et al., 2003; Gillett and Fyfe, 2013), which suggests a decline in the age-1+ population in both regions.

We demonstrated that decadal climate variation and environmental drivers significantly impact the recruitment of Arctic cod early life stages into the age-1+ population. The NAO and AO are forecasted to become more positive throughout all seasons (Gillett et al., 2003; Gillett and Fyfe, 2013), which our observations suggest could lead to an eventual decline in the age-1+ Arctic cod biomass in both the Beaufort Sea and Baffin Bay, the results of which would cascade throughout entire Arctic marine food webs (Pedro et al., 2023). However, we also show that the response of Arctic cod populations to environmental and climatological changes varies spatially across the Canadian Arctic and that 13%–42% of the total variation in age-1+ Arctic cod biomass could not be explained by the variables tested here, with less variation explained in Baffin Bay. Forage fish recruitment rates are typically highly variable and other mechanisms not tested here could be at play (e.g., current velocity, freshwater input, larval transport and retention, predation), especially at finer scales. In addition, examining these data on a finer spatial scale in relation to more localized environmental data may help elucidate more defined patterns for each area, and movement studies may provide a better understanding of how age-0 fish are influenced by surface currents and climate indices. Yet, we have shown clearly that variation in environmental drivers during the early life stages of Arctic cod affects the recruitment into the age-1+ population. The low pelagic fish biodiversity of the Canadian Arctic, coupled with a food web that is largely controlled by bottom-up processes, means that this region may be less resilient to changes in Arctic cod dynamics than more diverse regions. These changes would have both ecological and sociological impacts, because Greenland halibut (Reinhardtius hippoglossoides) and Arctic char (Salvelinus alpinus), the main fish species harvested in the Canadian Arctic, prey on Arctic cod (Vollen et al., 2004; Peklova et al., 2012; Spares et al., 2012). Shifts in Arctic cod dynamics would thus greatly impact both subsistence and commercial fisheries.

The acoustic-trawl datasets collected by the Canadian research icebreaker CCGS Amundsen are made available by the Amundsen Science program, which is supported through the Université Laval by the Canada Foundation for Innovation. These data from 2006–2016 are available on https://www.polardata.ca (CCIN Reference No. 12841) or the entire data set can be obtained by emailing [email protected]. The acoustic-trawl datasets collected by the FV Frosti are made available by the Department of Fisheries and Oceans Canada. These data can be obtained by emailing [email protected].

The supplemental files for this article can be found as follows:

Figures S1–S10. Tables S1–S2. DOCX

The authors thank C. Aubry for his help with field sampling and database management, W. Hiscock for his help with acoustic data cleaning, the technicians and students who participated in sampling, and all the crew members and technicians aboard the CCGS Amundsen and the FV Frosti for their invaluable help over the years. They are grateful to the Inuvialuit Game Council for supporting data and sample collection aboard the FV Frosti and for providing valuable input regarding sampling and logistics. In memoriam of Louis Fortier (1953–2020), a major contributor to the previous studies that inspired this one and to the provision of a large proportion of this data set.

Funding was provided by the Natural Sciences and Engineering Research Council (NSERC) through the Discovery Grant program, the ArcticNet a Network of Centres of Excellence Canada, the Beaufort Regional Environmental Assessment (BREA) (Aboriginal Affairs and Northern Development Canada), the Environmental Studies Research Fund (Natural Resources Canada, NRCAN), the Program of Energy Research and Development (NRCAN), the International Governance Strategy and other internal funding from Fisheries and Oceans Canada (DFO), and Crown-Indigenous Relations and Northern Affairs Canada (CIR-NAC). JH, JF, and MG were funded through NSERC and ArcticNet. MG was funded through CIR-NAC. AN, AM, and SG were funded through BREA, NRCAN programs, and DFO. Some of the data presented herein were collected by the Canadian research icebreaker CCGS Amundsen and made available by the Amundsen Science program, which was supported by the Canada Foundation for Innovation and Natural Sciences and Engineering Research Council of Canada. Research permits for work aboard the Amundsen were obtained through ArcticNet (https://amundsenscience.com/expeditions/past-expeditions/).

The authors have no competing interests to declare.

Conception and design: JH, JF, CB, AN, MG.

Acquisition of data: MLB, AM.

Contributed to analysis and interpretation of data: JH, JF, MG.

Drafted and/or revised the article: JH, JF, CB, AN, MLB, AM, SG, MG.

Approved the submitted version for publication: JH, JF, CB, AN, MLB, AM, SG, MG.

How to cite this article: Herbig, J, Fisher, J, Bouchard, C, Niemi, A, LeBlanc, M, Majewski, A, Gauthier, S, Geoffroy, M. 2023. Climate and juvenile recruitment as drivers of Arctic cod (Boreogadus saida) dynamics in two Canadian Arctic seas. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2023.00033

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: Four Decades of Arctic Climate Change: A Tribute to Louis Fortier