Arctic cod (Boreogadus saida) hatching in the Hudson Bay system: Testing of the freshwater winter refuge hypothesis

Arctic cod (Boreogadus saida) spawn in late fall and winter when the buoyant eggs rise to the surface or to the ice-water interface (Ponomarenko, 2000). These eggs hatch from January to early July depending on the Arctic or sub-Arctic region (Bouchard and Fortier, 2008, 2011). The eggs incubate for a minimum of 31–35 days, but the colder the temperature, the longer the incubation time (Laurel et al., 2018). The hatch length of Arctic cod larvae is usually not more than 6 mm (Aronovich et al., 1975). The larvae start to absorb their yolk sac and very soon after hatching, before the yolk is fully absorbed, start feeding on nauplii and small mesozooplankton in the surface layer (Andersen et al., 1994). They start swimming faster and gradually start hunting for larger prey when they are about 30–35 mm in length (Baranenkova et al., 1966; Bouchard and Fortier, 2020; Copeman et al., 2020).

Individuals who hatch early at low temperatures under the ice face relatively poor feeding and survival conditions. Before the spring bloom, food is sparce, and foraging is more difficult due to low temperatures and low light. However, early hatchers would benefit from a longer growth period. The early hatchers that do survive would reach a larger size in autumn than late hatchers (Fortier et al., 2006), which presumably provides a survival advantage during the following winter. Larger juvenile fish are better at surviving the low temperatures and scarce food during the winter months because of their larger lipid reserves (Copeman et al., 2020). Additionally, larger juveniles tend to migrate toward greater depths in the fall where they are less conspicuous to visual predators (Geoffroy et al., 2016). Winter hatching has been observed in Arctic seas strongly influenced by freshwater inputs, such as the Laptev Sea, the Beaufort Sea, and Hudson Bay (Bouchard and Fortier, 2011). While winter temperatures under the ice in marine waters are always at freezing point (–1.8°C for a salinity of 33), temperatures in waters influenced by under-ice river plumes are slightly warmer (for example –1°C for a salinity of 20). This difference has physiological implications on the early life dynamics of Arctic cod because 0°C–4°C is the optimal temperature for hatching success, growth rate and survival (Laurel et al., 2017, 2018; Koenker et al., 2018). The notion that some Arctic cod eggs drift to and hatch in brackish water leading to enhanced larval survival is called the freshwater winter refuge hypothesis (Bouchard and Fortier, 2011; Bouchard et al., 2015).

The large amount of freshwater entering the Hudson Bay means that potential thermal refuges under the ice in winter are plentiful (Déry et al., 2005). The geographic location of Hudson Bay at the southern margin of the Arctic Ocean system grants the bay features that match near-future projections for most of the inner Arctic Ocean. First, sea ice cover will become truly seasonal (no more multiyear ice; Serreze, 2011), and second, freshwater input should increase with further climate warming (Déry et al., 2016). These parallels make the study of the Hudson Bay marine ecosystem particularly useful, as observations made in the bay now might inform on the prevailing conditions in higher latitude ecosystems in the near future (Macdonald and Kuzyk, 2011).

Due to the elusive nature of the drift of Arctic cod larvae during their early lives (especially in ice-covered areas), understanding the functional relationships between their hatch location, their survival and the recruitment of juveniles remains challenging. In this study, we used two independent and complementary approaches to estimate the environmental conditions at the time and location of hatching. First, because both the eggs and early larval stages of Arctic cod are buoyant (Spencer et al., 2020), we used a high-resolution coupled biophysical model to simulate the larval trajectories from the catch location back to potential hatch locations. Second, we analyzed the microchemistry in the otolith core to estimate the salinity at the actual hatch location.

The Hudson Bay system comprises the Hudson Bay itself and the adjacent Foxe Basin, Hudson Strait, Ungava Bay and James Bay (Figure 1). It receives freshwater from over 42 rivers draining a large watershed extending over most of Canada and parts of the United States (Déry et al., 2005). The Hudson Bay is a shallow inland sea (maximum depth in the center of approximately 250 m), and the deepest part of the system is in Hudson Strait (maximum depth of approximately 1,000 m) which opens to the Labrador Sea in the North Atlantic. The Hudson Bay system is also highly stratified. The surface layer is formed by fresher and warmer water coming from the rivers, whereas the bottom layer is composed of denser marine water masses entering Hudson Bay from the Arctic via the Roes Welcome Sound and from the Atlantic Ocean along the northern coast of Hudson Strait (Drinkwater, 1988; Stewart and Lockhart, 2005; Defossez et al., 2010; Ridenour et al., 2021). The surface circulation is largely cyclonic and surface water is exported out of the Hudson Bay system via Hudson Strait (Ridenour et al., 2019b). Partly due to the large input of freshwater increasing the freezing point of its surface waters slightly, the Hudson Bay is seasonally covered with sea ice despite its low latitude (Prinsenberg, 1986). The arctic climate and the westerly winds prevailing over the region lead to the formation of the Hudson Bay ice cover, usually starting in November and lasting until late June or early July (Wang et al., 2012).

Five expeditions of the research icebreaker CCGS Amundsen were conducted in the Hudson Bay system between 2005 and 2018 (Table 1). Conductivity-temperature-depth (CTD) profiles performed with a Seabird SBE19plus® sensor were used to estimate surface salinities (average over 0–20 m) and separate stations between higher freshwater influence (HFW) and stations with lower freshwater influence (LFW), using the median of all stations and years as a threshold value.

All Arctic cod larvae captured during the above-mentioned expeditions, a total of 160, were included in this study. The larvae were caught with five types of obliquely trawled samplers: (1) a double square net, (2) a rectangular midwater trawl, (3) a multilayer sampler, (4) a large pelagic trawl, and (5) a hand-operated ring net. Technical specifications and deployment methods for the first four samplers are described in Bouchard and Fortier (2011). The fifth sampler was a 65-cm diameter aperture, 3.25-m long, 500-µm mesh net deployed from a zodiac or a barge in locations inaccessible to the ship near landfast ice or in estuaries. The ring net was deployed to a depth of around 10 m. All samplers were equipped with rigid, filtering cod-ends with a mesh size corresponding to that of the net.

All age-0 fish collected (larvae and juveniles less than a year old) were identified to the lowest taxonomic level possible on board. The early life stages of Arctic cod (Boreogadus saida) and ice cod (Arctogadus glacialis) are very similar morphologically and only distinguishable using genetics, otolithometry or size-at-date frequency distribution (Madsen et al., 2009; Bouchard and Fortier, 2011; Bouchard et al., 2013; Bouchard et al., 2016). Ice cod is a High Arctic species, and none of the 38 larval Gadidae collected in Hudson Bay in 2005 that had been analyzed genetically belonged to this species (Bouchard and Fortier, 2011). Therefore, in this study, all Gadidae were assumed to be Arctic cod. All age-0 Arctic cod (hereafter referred to as larvae), except two highly degraded individuals, were measured on board (standard length, SL) and preserved individually in 95% ethanol.

To include a measure of the time passed since ice break-up in the analysis, the variable “days of open water” (DOW) was included in the analysis. This variable was defined as the number of days from when total ice concentration, extracted from regional ice charts produced on the first Monday of each week by the Canadian Ice Service, dropped below 50% up to the catch date. The DOW was calculated at each individual station.

Several trace elements dissolved in the water where the larvae develop are incorporated daily into the carbonate matrix of the growing otolith and represent natural tags of what is found in the water (Campana, 1999). The concentration of strontium (Sr) in the otolith is known to be related to the salinity of the environment, as this element is found in lower concentrations in fresh and brackish waters of riverine origin and from sea ice meltwater (McMullin et al., 2017). Salinity at the hatch location can therefore be inferred from the concentration of Sr in the otolith core (Elsdon et al., 2008). Otoliths were extracted first for age reading and then for microchemistry analysis. To save time while conserving a representative sample, otoliths were extracted from a maximum of 10 larvae in each sampling year and 5-mm length class, and from 75% of the larvae in larger sub-groups. The lapilli were extracted from the preserved larvae under a stereomicroscope with polarized filter and mounted on microscope slides on their medial side using Crystal Bond® thermoplastic glue. Mounted left lapillar otoliths were observed under an Olympus BX51 light microscope at 1000× magnification, and polished on the medial side using 5 µm aluminum grit paper until all of the increments became differentiated at 1000× magnification. Images of the otoliths were taken using an Olympus DP70 digital camera coupled to a light microscope, and the number of increments (day rings) were counted to obtain the larval age. Other measurements of increment widths, core radius, and otolith radius were also taken using Image-Pro Plus® image analysis software.

For the 133 larvae aged by otolithometry, standard length and age correlated linearly (age = –7.26 + 4.251 SL, R² = 0.96, P < 0.0001, F_(1,130) = 3,201), hence the age of the remaining 25 larvae was estimated from their standard length by randomly selecting a value from the normal 95% confidence interval around individual fitted values predicted for each standard length. Then, in order to compare the standard length between all individuals that hatched at different times during the season, the standard length on one specific date was estimated. October 2 was selected as it was the latest day in our dataset when Arctic cod larvae were sampled. The standard length on October 2 (SL_2O) was predicted from the age that each larva would have been on that date (October 2 minus hatch date) using the aforementioned method of regression.

The open source, individual-based modeling framework Ichthyop uses physical and biological dynamics to simulate the distribution of fish eggs and larvae through time (Lett et al., 2008). The objective of the Lagrangian simulations was to estimate the salinity that the Arctic cod larvae may have experienced in a potential hatch location area. Ichthyop was run backward in time from the catch date to the hatch date as estimated from otolith ageing for 158 larvae. In each simulation, 1,000 particles were released at the location and time of the sampling to ensure that the variability in the estimation of the hatch locations resulting from the stochastic process of Lagrangian transport was captured. The time step used was 1,080 s (18 min). The salinity and temperature values of the surface layer of all of the potential hatch locations were recorded, and the median salinity is considered the predicted salinity hereafter.

The parameters used for the Ichtyop simulations are shown in Table 2. The three-dimensional advection forcing fields, as well as the temperature and salinity information, came from simulations of the Nucleus for European Modelling of the Ocean version 3.4 (Madec and NEMO System Team) set up with the Arctic and Northern Hemisphere Atlantic 1/12° resolution configuration (ANHA12; Hu et al., 2018; Ridenour et al., 2019a, 2019b). The ice module used was the Louvain-la-Neuve Ice Model version 2 (LIM2; Hunke and Dukowicz, 1997). The forcing outputs are provided at 5-day averages. There are 50 vertical levels with the surface layer having a thickness of 0.5 m, then increasing with depth. The atmospheric forcing applied to the ANHA domain is the Canadian Meteorological Centre’s global deterministic prediction system reforecasts (Smith et al., 2014) from 2002 to 2018. The model was initiated with 3D fields of temperature, salinity and horizontal velocities and 2D fields of sea-surface height and sea ice concentrations taken from the GLobal Ocean ReanalYsis and Simulations (GLORYS2v3) provided by Mercator Ocean (Masina et al., 2017). River input was provided by the hydrological forcing found in Dai and Trenberth (2002) and Dai et al. (2009), whereas solid and liquid glacial discharges for Greenland were provided by Bamber et al. (2012). Note that river runoff forcing was repeated from 2007 to 2018, as 2007 was the most recent year for which such data were available (see Ridenour et al., 2019a, for more details).

Horizontal diffusivity is not used in backtracking simulations because diffusive processes are stochastic and not reversible in time. But advection fields from high-resolution regional circulation models such as ours (see below; Hu and Myers, 2013) have been shown to capture the mean pathways of particles during backward tracking simulations adequately (Gillard et al., 2016; van Sebille et al., 2018). Previous studies have used backtracking with Ichthyop to assess the environmental conditions experienced by Sargassum algae floats in the North and equatorial Atlantic (Putman et al., 2018; Michotey et al., 2020) and for studies on the Mediterranean coastal fish hatch distributions (Calò et al., 2018).

To add an independent estimate of the salinity at the hatch location, we measured Sr concentration in the core of the left lapillus. A low Sr concentration in the core indicates that the salinity at the hatch location was low, and vice versa. The otoliths from the Arctic cod captured in 2005 were analyzed for otolith microchemistry, as reported in Bouchard et al. (2015), and destroyed in the process. The Sr concentration in the otolith cores of the 2005 individuals could not be included in this study because different ablation methods were used and the Sr concentrations are not directly comparable. Therefore, the otolith microchemistry of 87 individuals were analyzed from the 2010, 2017 and 2018 samplings. The left lapillus was analyzed but to verify the reliability of the results, the left and the right lapillus were analyzed for nine of the individuals for a total of 96 otoliths analyzed over two sessions at the Laboratoire des matériaux terrestres at Université du Québec à Chicoutimi, Chicoutimi, Canada.

The otoliths had a small average radius of 11.2 µm, which presented a risk of losing them during preparation; therefore, they were handled minimally. The otoliths were transferred from their individual slide used for the 1000× magnification light microscope to another slide to be used in the mass spectrometer. Crystal Bond® thermoplastic glue was again used to bind the otoliths to the slide. The polishing done during age-reading was sufficient, and no further polishing was needed to prepare the otoliths for laser ablation.

The Laser Ablation Inductively Coupled Plasma Mass Spectrometry (LA-ICP-MS) analyses were undertaken using a Resonetic Excimer 193 nm ArF laser coupled to an Agilent model 7900 ICP-MS (Agilent, Mississauga, Ontario, Canada). Using image analysis software, the laser was positioned at the otolith edge so that a transect from one edge to the other, passing through the core, was made (Figure S1). For all analyses, the laser beam diameter was set to 15 µm, moving at a speed of 2 µm s^–1, at a frequency of 20 Hz and a laser fluence of 5 J cm^–2. Each analysis was normalized to an expected calcium content (⁴⁴Ca) in otolith matrix (38.02% wt; Campana, 1999). A NIST-610 glass reference material was measured regularly for 60 s after every eight samples to estimate element concentrations and correct for any analytical drift. Additional reference materials USGS, MACS-3, and GP-4 were also measured for quality control. Carrier gas background without ablation was measured for 30 s before analysis and 15 s after analysis and then subtracted from the observed values.

The transect through the otolith was recorded as a series of concentrations with a corresponding time stamp. To establish the correspondence between sections of the transect with regions in the otolith, the time when the laser passed through the center of the core was recorded. This time and the size of the core measured during age determination were used to determine the part of the transect that had Sr measurements belonging to the core. Measurements of the lapillus core showed that it had an average radius of 8.2 µm. As the laser travelled at 2 µm s^–1, the concentration of the core was the average of the concentrations at 2.05 s on either side of the center of the core.

Concentrations of Sr from the ICP-MS signal were estimated using Iolite (Paton et al., 2011), a free add-on to Igor Pro software (Wavemetrics Inc, Portland, OR, USA). Data integration was performed using the trace element IS procedure in Iolite. To assess the validity of observed concentrations, the limit of detection of the element was calculated as three times the standard deviation of the gas blank divided by the sensitivity of the signal (Lazartigues et al., 2014).

We tested the relationships between the response variables hatch date (HD) and standard length on October 2 (SL_2O) with the major environmental variables experienced by Arctic cod larvae using linear mixed-effects models. The explanatory variables tested were the catch date (CD), the days of open water (DOW) and two distinct indices of the salinity at hatch: the predicted salinity (PS) at the simulated hatch time and location, and the Sr concentration in the otolith core (Sr). However, in order to keep both estimations independent and complementary, we performed the linear mixed model analyses separately for each index (Table 3). These models also included all of the possible interaction terms.

Because sampling was carried out in different years with different environmental conditions and at a different time of year, the year (Y) variable was included in the model as a random effect. Interaction terms between the explanatory variables and random slopes were also tested. Explanatory variables and interactions were removed sequentially from the linear mixed-effects model and left out if they did not have a significant effect on the slope.

Prior to fitting the models, all of the variables were tested for normality and transformed accordingly. Neither HD nor SL_2O were distributed normally, hence HD needed a quadratic transformation (HD²), while a log transformation was needed for SL_2O (log SL_2O). Moreover, residuals were also tested after each test to verify that the assumptions of normality and homoscedasticity were satisfied.

A generalized linear model was used to test the regression between the Sr concentration in the otolith core and the PS. All analysis was done on R 4.0.4 with the "lmerTest" package for the mixed-effects modeling.

Most Arctic cod larvae hatched in May and June, although 8% of them hatched in January and February, some as early as January 9. Surface salinity (0–20 m) at the catch location ranged between 17.5 and 32.8. Stations with surface salinity below the median (30) were considered to have a high freshwater influence while those with surface salinity >30 were considered low freshwater influence. About two thirds of HFW stations were located in the southern Hudson Bay and near river mouths (Figure 2a). About 70% of the early hatchers (hatch date between January and April) were captured at HFW stations with surface salinities below the second quartile (Figure 2b). SL_2O was significantly higher at HFW stations (36.3 mm) than at LFW stations (32.2 mm; t(104) = –3.0, P = 0.004; Figure 2c). At the stations with the lowest surface salinity, the larvae had some of the largest SL_2O, but there were also several individuals with very low SL_2O. Therefore, the larvae in this group went against the trend of decreasing SL_2O with increasing surface salinity. The results showed some correlation between salinity and SL_2O, while highlighting the need to bridge the gap between environmental conditions at the catch and hatch locations.

In 2005 and 2010, sampling was carried out later in the season. The larvae caught were older, and therefore the simulated backward trajectories were longer. We found rather large direct linear distances of 182.0 km in 2005 and 225.7 km in 2010 between the catch locations and the simulated hatch locations (Figure 3a–f). During Lagrangian advection simulations, even particles starting close by each other can be transported by different prevailing currents. Therefore, the longer tracks in 2005 and 2010 resulted in several potential hatch areas (Figure 3a and d). As a result, the predicted surface salinity among the potential hatch areas varied by up to 5 salinity units, especially when they were close to a river plume where large variations in salinity occurred within a small distance. On the contrary, the estimated hatch locations for most larvae in 2017 and 2018 were much closer to the catch locations, on average 52.1 and 59.9 km, respectively (Figure 3g–l). The potential hatch locations determined by the backward simulations were all similar for the 1,000 particles released on each run (Figure 3g and j).

While particles tracked by the model could die off if they entered water with temperatures above (or below) some lethal thresholds, this die-off only happened for some of the larvae coming from one station in 2005 when temperatures exceeded 11°C along the southeast coast of Hudson Bay. Temperatures at the hatch locations varied only slightly, within a single degree Celsius, as winter temperatures are quite stable. Nonetheless, surface temperatures at the hatch locations were the lowest in 2018 and the highest in 2010. In 2005, the temperatures at the hatch locations were similar to those in 2017 and 2018, even though the median hatch date was 2 months earlier (Table 4).

Larvae that hatched earlier tended to experience lower predicted salinity at their hatch locations (DF = 3, ΔAIC = 38.4, P < 0.001; Figure 4a). However, the mixed-effects model showed that this relationship varied substantially among years, being totally reversed in 2018 and almost absent in 2017 (Figure 4d). Year had a very significant random effect (DF = 2, ΔAIC = 30, P < 0.001; Figure 4d). The time of sampling affected the slope of the relationship between hatch date and predicted salinity. Catch date and the days of open water were significant variables in the mixed-effects model (DF = 2, ΔAIC = 15.2, P < 0.001 for CD; Figure 4b, and DF = 2, ΔAIC = 15, P < 0.001 for DOW; Figure 4c). Early hatchers were more likely to be captured when sampling was done later in the year after the sea ice had melted. The interaction between catch date and days of open water was also significant (DF = 1, ΔAIC = 16.9, P < 0.001). When sampling was carried out in September and October, more time had passed since the 50% of ice break-up, and this interaction term explained a small but significant amount of variation in the mixed-effects model.

Larvae that had a larger SL_2O had also hatched in lower predicted salinity (DF = 3, ΔAIC = 40.87, P < 0.001; Figure 5a). The catch date was a significant fixed effect in explaining variation in SL_2O (DF = 2, ΔAIC = 11.47, P < 0.001; Figure 5b), as was the interaction term between catch date and days of open water (DF = 1, ΔAIC = 13.46, P < 0.001). Larvae that were projected to achieve a larger size by the end of summer were often captured at stations sampled later in the year when the sea had been ice-free for more than a month (DF = 2, ΔAIC = 11.52, P < 0.001; Figure 5c). Again, there was a strongly significant random effect from the year variable, with 2018 showing opposite slopes from 2005, 2010, and 2017 (DF = 2, ΔAIC = 29.61, P < 0.001; Figure 5d).

The Sr concentration in the otolith cores ranged from 1,373 to 4,764 ppm (Figure 6). We used a concentration of 2,000 ppm or lower in an otolith core to indicate that the fish spent its first days in brackish water (approximately salinity below 30) (Morissette et al., 2016). A few otoliths had Sr concentration of less than 2,000 ppm, indicating that a handful of larvae hatched in waters affected by freshwater plumes. Only eight otolith cores had a Sr concentration above 3,000 ppm, indicating that hatching in more saline waters was rare. No significant relationship was found between Sr concentration in the otolith core and the estimated salinity at the hatch location (p = 0.196; Figure 6), neither for the pooled data nor for separate years. The most common concentrations of Sr in the otolith cores were between 2,500 and 3,000 ppm, suggesting that most fish hatched in a salinity range of 27 to 32. The predicted salinities from the modeled hatch locations were in this range. However, there was also a group of otolith cores with Sr concentrations below 2,000 ppm and no larvae with a low predicted salinity at the hatch location. About 30% of the larvae tested showed mismatches between predicted salinity at the hatch location and the Sr concentration in their otolith cores (Figure 7). However, close inspection of these results reveals that individuals with a high predicted salinity but low Sr were transported very short distances and were mostly from the same area in the northwest part of Hudson Bay. In many cases, areas of low salinity were close to the larvae tracks (Figure 7a). Individuals with high Sr and low predicted salinity were located along the eastern coast of the bay, and at one station along the southern shore of Hudson Strait (Figure 7b). Individuals showing consistent metrics could be found anywhere else (Figure 7c). Notably, the individuals with the largest mismatches between modeling predictions and Sr concentration results were all captured in 2018.

Larvae with low Sr concentrations in the otolith core generally hatched earlier (DF = 2, ΔAIC = 2.3, P < 0.005; Figure 8a) than those with a high Sr concentration, but the variance explained by Sr concentration was small. The catch date had a significant positive relationship with the hatch date (DF = 1, ΔAIC = 38.7, P < 0.001; Figure 8b). The days of open water was also an important explanatory variable for the variation in hatch date (DF = 2, ΔAIC = 2.4, P < 0.05; Figure 8c). When larvae were captured soon after ice break-up, the probability of catching late hatchers was higher. The interaction term between Sr concentration and days of open water (DF = 1, ΔAIC = 2.5, P < 0.05) also showed that the probability of catching a larva with a higher Sr concentration increased as DOW decreased, meaning that the larvae sampled were more likely to have hatched in higher salinity. The random effect of the year of sampling was also significant (DF = 1, ΔAIC = 75, P < 0.001; Figure 8d) with the earliest hatching larvae being caught in 2010 and the latest hatching ones in 2018.

Finally, no significant relationship could be found between Sr concentration in the otolith core and SL_2O by the linear mixed-effects model (DF = 1, ΔAIC = 1.29, p = 0.40). The catch date was the only significant fixed effect in explaining the variation in SL_2O (DF = 1, ΔAIC = 45.83, P < 0.001), and year as a random effect was also significant (DF = 1, ΔAIC = 47.13, P < 0.001).

Arctic cod in the Hudson Bay system tend to hatch earlier than in other regions of the Arctic (Bouchard and Fortier, 2011), a feature that has been attributed to the influence of the large volume of freshwater that enters the bay from summer to early winter (Déry et al., 2011). By analyzing data at the scale of the Hudson Bay system, we were able to link early hatching larvae that survived until sampling with potential hatch locations influenced by freshwater input. We made this link first by simulating the path of individual larvae to their potential hatch locations and second by analyzing the otolith microchemistry for indices of freshwater influence. The Arctic cod larvae that hatched earlier were almost always traced back to surface salinity below the median and had lower Sr concentrations in their otolith core. Most of these early hatching larvae were found to have hatched in the southeastern part of the bay where the surface salinity was lower.

The backtracking experiments run from the catch date to the hatch date of individual larvae show that many trajectories end up in coastal areas, not far from river estuaries, especially along the east coast of Hudson Bay and the south coast of Hudson Strait. Another important outcome of our modeling experiments is that in Hudson Bay larvae are not necessarily advected over great distances, especially those that hatched after April. In 2010, for example, the larvae captured in the vicinity of the Nelson River estuary were traced back to a hatch location further into the estuary, meaning that the larvae must have drifted against the large scale cyclonic coastal current that persists in this area throughout the open water season (Ridenour et al., 2019b). This movement against prevailing currents shows the influence that local currents can have on trajectories for Arctic cod larvae and highlight the importance of using forcing fields from a regional circulation model with a spatial resolution at least as high as ours (1/12th of a degree or approximately 3.5–5.5 km). The finding that larvae remain confined to the vicinity of their hatch location also warns that they are vulnerable to sudden changes in local environmental conditions (e.g., temperature) due to, for example, river runoff regulation.

The intrinsically chaotic nature of Lagrangian experiments means that small divergences in the modeled trajectories accumulate over time (Christensen et al., 2007). This accumulation essentially affected the hatch location of the larvae that had the longest trajectories, especially those reaching the northern part of Hudson Bay where currents from the Foxe Basin, Hudson Bay and Hudson Strait meet (Figure 1). For example, simulated trajectories from 2005 showed that hatch location estimates for particles released at the same location could end up hundreds of kilometers apart. Any calculation for estimating the likelihood of one hatch location over another would still carry the same model biases. However, we also provide a complementary index for the salinity at the hatch location: Sr concentration in the otolith cores.

Given the assumptions and uncertainties inherent to both Lagrangian modeling and the use of Sr concentration as an index of freshwater exposure during development, finding a direct correlation between the predicted salinities at the estimated hatch locations and Sr concentration in the otolith core of the sampled Arctic cod larvae was unlikely. As expected according to otolith microchemistry, more larvae should have hatched in a low salinity environment than predicted by the Lagrangian simulations. Noticeably, most discrepancies between the simulations and the microchemistry results were observed in 2018, a year during which precipitation was less than usual over the Hudson Bay System (Lukovich et al., 2021). The regional circulation model was highly resolved vertically (layers of about 0.5 m close to the surface), sufficiently to simulate surface plumes from either riverine or sea ice melt origin. However, the river runoff forcing fields used in 2018 were repeated from 2013 (see details in Ridenour et al., 2019a). Hence, the simulated salinity fields were likely less accurate in 2018 at the small scale relevant for our study of larval transport.

Nonetheless, the Lagrangian experiments and the Sr analyses provide coherent insight into the fate of Arctic cod larvae. The linear mixed models separately regressing the hatch date to the predicted salinity and the Sr core concentration predicted a significant positive interaction with both explanatory variables. These statistical analyses show that the larvae traced back within the lowest salinities (<27) and the ones whose otoliths presented the lowest Sr concentration (<2,000 ppm) both hatched the earliest. Conversely, the larvae traced back within the highest salinities (>31) and the ones whose otoliths presented the highest Sr concentration (>2,800 ppm) both hatched later. In the former case, the larvae reached the largest standard length, whereas in the latter case, they formed the group with the smallest standard length, although the correlation between Sr concentration in the otolith core and SL_2O was not statistically significant overall. A smaller standard length is assumed to be a disadvantage for overwintering, as smaller larvae have lower fat reserves (Copeman et al., 2020) and are more likely to be preyed upon or to die of starvation. For salinities and Sr concentrations in the middle of their respective ranges, there was a mix of early and late hatchers. These independent results paint a coherent picture: Arctic cod eggs in lower salinity will be more likely to hatch early, whereas eggs in high salinity are more likely to hatch later in the season and have potentially lower chances of survival during the following winter.

At the heart of the freshwater winter refuge hypothesis lies this simple relationship between salinity and temperature: During winter, saltier waters can withstand colder temperatures (before freezing), thus offering a lower quality habitat for Arctic cod larvae than fresher water masses. In the Greenland Sea, sea surface temperatures between –1.7°C and –1.0°C were associated with poor larval survival, whereas temperatures slightly above 0°C were associated with good survival (Fortier et al., 2006). In the laboratory, the survival of late-stage larvae was better at 2°C than at 0°C; as their metabolic rate and cardiac activity increased with increasing temperatures, Arctic cod larvae could increase both their lipid reserves and mass-to-length ratio (Koenker et al. 2018).

In our data set, the salinities at the predicted hatch locations ranged between 25.1 and 32.7. This difference in salinity could result in a temperature 0.5°C warmer within a thermal refuge. Assuming a Q₁₀ coefficient of 2 for the Arctic cod larvae metabolism and growth (i.e., rates double with a change of 10°C; Laurel et al., 2017), an increase of 0.5°C gives a rate (R) change of $R2=R1×20.510=R1×1.03$⁠, that is, a 3% increase of the daily growth rate. As small as that increase may seem, in only 23 days a larva will be twice as large as another that would not have benefited from this slightly higher growth rate. Hence, even a modest temperature increase of 0.5°C could result in a significant advantage for early hatchers that have on average about an extra month to grow. Early hatchers will be in better condition and will increase their chances of surviving the following winter to recruit into the age-1 population the next year.

The relative importance of developing in warmer waters and hatching early can be seen in our results. The predicted temperature at the estimated hatch locations was the highest in 2010, and among 2010, 2017, and 2018, the 2010 larvae achieved the largest standard length. However, the SL_2O of the 2005 larvae were the largest by far, even though the estimated temperatures at hatch were similar to those of 2017 and 2018. This apparent discrepancy is likely due to the hatching being earlier by 1–2 months for the larvae sampled in 2005. Our results show that while even a slightly warmer ambient temperature enhances larval growth rate, in the long run having several more weeks of growth by hatching early is a surer way to gain a larger size before winter.

In this particular Arctic ecosystem, warm temperatures can also have a detrimental impact on Arctic cod larvae. In Hudson Bay, hatching early might indeed provide larvae the time to grow enough to avoid lethally high temperatures later in the season (Bouchard et al., 2020). The surface temperature in this highly stratified system often rises above 8°C due to the spring freshet bringing warmer river water right when the thermal stratification takes place in summer (Ridenour et al., 2019a). As Arctic cod larvae are most sensitive to heat and salinity changes during their first 6 weeks posthatch (Spencer et al., 2020), early hatching in Hudson Bay could improve their chances of survival if they have already passed this critical 6-week period when the surface waters become warmer. This possibility was partly confirmed by the high simulated mortality reached at one location in 2005 where the upper temperature threshold was exceeded.

While ensuring mostly beneficial effects from their ambient temperature, surviving early hatchers will also benefit from prime access to their prey. Even when hatching months ahead of the short-lived spring bloom in May–June (Barbedo et al., 2020), Arctic cod larvae can rely on scarce yet suitable zooplankton prey, essentially composed of small-sized copepod species (e.g., Pseudocalanus spp., Microcalanus spp., Oithona spp.) or nauplii and copepodite stages of larger, lipid-rich copepod species (e.g., Calanus spp.). Recent studies have shown that some primary production occurs as soon as February in Arctic seas, even below the ice cover (Randelhoff et al., 2020). This early production implies that a small proportion of the lipid-rich overwintering copepods might leave diapause during winter to join the smaller ones that remain active (i.e., grow and reproduce) close to the surface, depending on environmental and physiological conditions (Hobbs et al., 2020). Hence, some small copepods (e.g. Pseudocalanus spp., Oithona spp.) and nauplii from larger species (Calanus spp.) would likely be available for early hatching larvae to feed on. The strong tides in Hudson Bay result in flaw leads in the ice forming semi-diurnally in coastal regions during winter. These leads are biologically active regions that allow some light to enter the water column (Stadnyk et al., 2019) allowing Arctic cod larvae to detect their prey (Bouchard and Fortier, 2008).

Our results show that small changes in salinity and temperature could result in significant changes in hatch date, summer survival, and the eventual prewinter condition of individual Arctic cod larvae. The current trends in the evolution of the Arctic sea ice cover include a shrinking in surface and thickness, later formation and earlier break-up. Hence, the variability in local conditions of salinity and temperature in the surface waters is expected to increase during the transitional seasons, especially spring (Wassmann and Reigstad, 2011). Consequently, an additional advantage for Arctic cod larvae hatching early in winter could be that they develop within a more predictable and consistent environment in the foreseeable future. The current level of variability we observed for this study can help to address this issue.

We sampled Arctic cod in different years, different months, and different stages in their larval development. While this sampling allowed us to cover a large range of environmental conditions and biological responses, it also showcased the high interannual variability inherent to Arctic marine ecosystems. In the spring of 2017 and 2018, many newly hatched larvae (some of them only a week old) were captured. In contrast, in 2005, sampling was carried out in late summer at the end of the open water period. By this time of the year, the larvae captured had already managed to survive the summer. The hypothesis that hatching early gives larvae a survival advantage implies that early hatchers would be more abundant in the nets deployed in 2005 than in 2017 and 2018. In 2017 and 2018, some of the larvae that were captured would not likely have survived to the end of the summer. In the Northeast Water (northwestern Greenland Sea), Michaud et al. (1996) found a much higher proportion of larger Arctic cod larvae (>8.5 mm) when sampling in late June to early August than when sampling earlier in late May through early June. The authors concluded that the higher water temperatures in the later months provided better survival conditions for the late hatchers (<0°C in May and June; 0°C –4°C in July and August). Assuming that the larvae had already passed their most vulnerable 6-week period, higher sea surface temperatures later in the season would enhance the growth rate and result in a higher standard length (Koenker et al., 2018). However, in the present study, we know that the larger larvae captured later after ice break-up had hatched earlier than the small larvae captured in early June. This finding is evidence that the timing of sampling might bias the larval developmental stages that are captured, with smaller larvae more likely to be caught earlier in the year. This possibility warrants further investigation and should be considered carefully when interpreting multiyear and multiseasonal data.

Another possible explanation for the relatively large amount of early hatchers captured in 2005 is that 2005 was an exceptional year for the amount of freshwater runoff that flowed into the Hudson Bay in both summer and winter (Ferland et al., 2011; Déry et al., 2011; Déry et al., 2016; Stadnyk et al., 2019). Precipitation was higher than the 30-year average in each month of 2005. In contrast, in 2018, precipitation was lower than average in 10 months of the year according to the HdryoGFD (Global Forcing Data for Hydrology; Lukovich et al., 2021). Although there were individuals that hatched even before the spring freshet, a larger number of Arctic cod larvae in 2005 hatched in March and April, just as freshwater started flowing into the Hudson Bay system. The inter-seasonal and interannual variability in Arctic cod larvae dynamics needs to be better understood and more closely monitored if we aim to predict the impacts of projected future environmental conditions on the population dynamics of Arctic cod populations.

Owing to the large freshwater runoff into the Hudson Bay system, stratification is high and the surface layer salinity remains relatively low even in areas quite far from any river outlets (Saucier et al., 2004). Additionally, areas far from the coast get freshwater input when all of the sea ice in the Hudson Bay melts within a few weeks in summer (Macdonald and Kuzyk, 2011). The Hudson Bay is the southernmost Arctic sea where sizable populations of Arctic cod thrive, and the environmental conditions described above have been rather unique until now. For example, to the west, the Chukchi Sea is a similarly shallow and seasonally ice-covered sea, but it does not receive a comparable inflow of freshwater. There, age-0 Arctic cod are found more commonly along the ice edge in northern regions with colder, more saline conditions and early hatching has not been observed (Kono et al., 2016; De Robertis et al., 2017). Along the Russian coast, the Laptev Sea is a shallow Arctic sea with seasonal sea ice cover and considerable freshwater inflow that is dominated by its largest river, the Lena. Contrary to the Chukchi Sea, Arctic cod larvae in the Laptev Sea have been observed to hatch as early as December of the previous year within the plume of the Lena river and to reach large pre-winter lengths (Bouchard and Fortier, 2008, 2011; Mishin et al., 2018).

More freshwater is expected to flow into the Arctic overall (Déry et al., 2016) and the ice-free period is expected to lengthen (Andrews et al., 2018). The current Hudson Bay ecosystem may provide a potential example of the destiny of Arctic ecosystems further north. An essential factor to consider, however, when assessing how under-ice freshwater dynamics might impact Arctic cod larvae survival in the near future is that environmental conditions may become more variable. Models have shown that there will be a higher volume and more regular input of freshwater into the Arctic through the year, owing to an earlier spring freshet and increase of water input throughout autumn and winter; however, this trend is projected to be stronger in Eurasia than in North America, and river regulation will dampen this trend (Stadnyk et al., 2021).

As our results have shown, the impacts of environmental change will differ from one location to another during a single year and within a single location amongst successive years. Moreover, the responses of Arctic cod population will vary at different life stages. Earlier ice break-up has been shown to give some advantage to Arctic cod larvae survival and subsequent recruitment into the age-1 population (Bouchard et al., 2017). After the precarious first few weeks, higher sea surface temperatures are a good predictor of feeding success in larvae of all sizes (Michaud et al., 1996). However, a long-term reduction in days of ice cover has also been linked to a reduction in the amount of larger sized Arctic cod consumed by predators, indirectly indicating a decrease in abundance (Gaston et al., 2012; Provencher et al., 2012; Chambellant et al., 2013). Although the effects of climate warming in the Arctic might be deleterious in the long term, predicting when and where a tipping point toward detrimental conditions for Arctic cod populations might occur remains a challenge.

In this study, we showed that among Arctic cod larvae that survived to be sampled, those that hatched early tended to do so in areas of low surface salinity and to achieve a larger size by the end of summer. The reverse was true for the late hatchers, and the picture was variable for individuals in the middle of the pack. We reached these conclusions by combining extensive direct measurements of the Arctic cod larvae and of environmental conditions at the sampling location, numerical simulations with a high-resolution regional circulation model and microchemistry analysis of the otoliths.

We also showed that the interannual variability in both environmental conditions and properties of Arctic cod larvae should be considered carefully when combining multiyear data sets. This observation is important, as data remain scarce in the Arctic, which leads to combining data of diverse origins. The effects of climate change make conditions in the Arctic more variable from year to year, which justifies continued sampling of ichthyoplankton and its biotic and abiotic environment, while increasing the research effort to extend the periods of the year covered by sampling. Wintertime under-ice processes as well as transitional periods (early spring, late fall) are increasingly recognized as priority targets for future Arctic marine research. In this study, we have shown their importance for the hatching and survival of Arctic cod larvae.

The data set within this article has been archived in the Canadian Watershed Information Network (CanWIN) with DOI: https://doi.org/10.34992/80NT-WC67.

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

Figure S1. Left lapillus of Arctic cod larvae.

Table S1. The gauged rivers in the Hudson Bay system, adapted from Stadnyk et al. (2021).

Professor Louis Fortier (Université Laval, Canada) was the PhD supervisor of SS and contributed to study conception and design and to data analysis and interpretation until his passing in October 2020. We express our gratitude for his contribution to this manuscript and mentorship offered to SS. Thanks also to the crew and scientists on the CCGS Amundsen, especially Cyril Aubry for help with field work. Cyril Aubry also helped with lab work as did Anne-Lise Fortin with otolith microchemistry lab work. We would like to thank David Babb for extracting data from regional ice charts and Andrew Tefs for providing information on freshwater input into the Hudson Bay. This work is a contribution to the scientific program of Québec Océan and of the CNRS & Université Laval Joint International Laboratory Takuvik (UMI3376).

This research is a contribution to the BaySys project funded by Natural Science and Engineering Council of Canada (NSERC) and Manitoba Hydro.

The authors have declared that no competing interests exist.

Contributed to conception and design: SS, ID, LF, FM, CB.

Contributed to acquisition of data: SS, ID, PGM, PS.

Contributed to analysis and interpretation of data: SS, LF, FM, CB, PS.

Drafted and/or revised the article: SS, FM, CB with input from all authors.

Approved the submitted version for publication: SS, ID, PGM, PS, FM, CB.

How to cite this article: Schembri, S, Deschepper, I, Myers, PG, Sirois, P, Fortier, L, Bouchard, C, Maps, F. 2021. Arctic cod (Boreogadus saida) hatching in the Hudson Bay system: Testing of the freshwater winter refuge hypothesis. Elementa: Science of the Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2021.00042

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

Associate Editor: Laurenz Thomsen, School of Engineering and Science, Department of Earth and Space Sciences, Jacobs University Bremen, Germany

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: The Hudson Bay System Study (BaySys)