Effects of early life history traits and warming on Arctic cod prewinter length and recruitment

The Arctic Ocean and its coastal seas are changing fundamentally with climate warming (AMAP, 2019; Steiner et al., 2021). Declining sea-ice cover, increasing surface temperatures, and changing ocean circulation are strongly impacting the functioning of Arctic marine systems (Hunt Jr et al., 2016). Sea-ice decline is of particular concern, as it may contribute disproportionally to large-scale loss and fragmentation of sea-ice habitats and induce shifts in ice-associated species phenology that might lead to recruitment failure and population decline (Søreide et al., 2010; Brandt et al., 2023).

The Arctic cod (Boreogadus saida) plays an important role in Arctic marine ecosystems due to its high biomass and its importance as a food source for seabirds, fish, and marine mammals (Mueter et al., 2020; Geoffroy et al., 2023). Arctic cod life cycle is intimately connected to sea ice, suggesting a particular vulnerability of this species to changes in sea ice, particularly for the early life stages (Geoffroy et al., 2023). The spawning season of Arctic cod is very long, extending from November to April (Geoffroy and Priou, 2020). Consequently, eggs are spawned and generally develop under the ice at sub-zero temperatures (Ponomarenko, 2000). A long spawning season ensures an increase in offspring survival under fluctuating environmental conditions by spreading the environmental risk among multiple spawning events (Hočevar et al., 2021; Hutchings, 2021). Within a region and a year, larvae originating from different spawning events will have different hatch dates, and thus experience different environmental conditions during their development, which in turn will influence their growth and survival (Bouchard and Fortier, 2008; 2011; Bouchard et al., 2021).

As the Arctic ice-free season becomes longer and the summer warmer, Arctic cod is expected to mature earlier, potentially leading to a shift in the timing of their spawning (Nahrgang et al., 2016). Early spawning may allow a longer growth period for winter and spring hatchers (Bouchard et al., 2017), but at the expense of increased predation exposure for the eggs, which develop slowly under the low temperatures of winter. Arctic cod eggs incubate for up to 79 days at –1.5°C and for only 35 days at 1.5°C (Aronovich et al., 1975; Laurel et al., 2018). While the egg stage can tolerate temperatures between –1.5°C and 3°C, the optimal embryonic development seems to occur just above 0°C, with a steep decline in hatch success observed above 2°C (Dahlke et al., 2018; Laurel et al., 2018). Arctic cod larvae can tolerate higher temperatures than the eggs, but their general fitness and survival are severely impaired above 5°C (Koenker et al., 2018a; Koenker et al., 2018b). While a moderate increase in temperature initially appears advantageous for Arctic cod by accelerating embryonic development and boosting larval growth, it carries the risk of narrowing the optimal thermal window for the development of early life stages.

The vulnerability of early life stages to climate warming and changing environmental conditions makes the first growing season critical, as larvae need to reach a level of fitness (e.g., length, weight, and lipid reserves) by the end of summer to ensure survival over their first winter (Copeman et al., 2022). Arctic cod larvae hatch at an average length of 6 mm (full range, 3.5–7 mm; Sakurai et al., 1998; Laurel et al., 2018), and by the end of summer the majority metamorphosize into juveniles, that is, adult-looking individuals with lengths between 25 mm and 35 mm and fully developed swimming fins (Ponomarenko, 2000; Walkusz et al., 2011). With increased swimming capacities, juvenile Arctic cod are migrating into deeper water by the end of summer (Geoffroy et al., 2016). However, an unknown proportion of juveniles, presumably the smaller individuals, remain in surface waters and become ice-associated in drifting sea ice (David et al., 2016). With significant changes in sea-ice cover and timing of ice break-up over the past decades, severe consequences are anticipated in the water masses such as increased temperature and changes in primary and secondary production that might affect the development of early life stages of Arctic cod (Mueter et al., 2020; Geoffroy et al., 2023). While increased primary production has been observed over the past decade in the Arctic Ocean, the future of Calanus copepods, the preferred prey for age-0 Arctic cod (Bouchard and Fortier, 2020), is less clear (Feng et al., 2018; Ershova et al., 2021; Brandt et al., 2023). A potential decline in calanoid copepods could result in low availability of a high-energetic food source for Arctic cod, for which they might not be able to compensate. Our understanding of this chain of consequences, however, is hampered by limited observations in the ice-covered Arctic Ocean, particularly in wintertime.

Our understanding of ecological processes under the ice-covered Arctic Ocean, which is rather poor, can be extended with the use of numerical models, allowing new scientific hypotheses to be tested with limited data at hand. Coupled ocean circulation models with sea-ice components have recently been used to simulate egg and larval Arctic cod dispersal (Huserbråten et al., 2019; Eriksen et al., 2020; Vestfals et al., 2021) and to infer spawning and hatching areas by backtracking larval dispersal from their sampling locations (Eriksen et al., 2020; Schembri et al., 2022). However, these coupled ocean circulation-transport models lacked important biological features of dispersing larvae, even though Vestfals et al. (2021) used a laboratory-fitted temperature-dependent regression in their model to estimate the growth of dispersing larvae. In addition to temperature, other factors are crucial for the growth dynamics of larvae, such as food availability. Using a bioenergetic model of larval Arctic cod growth, David et al. (2022) showed that the quantity of ingested food and the energetic content explained higher variability in modeled larval length than temperature alone. Moreover, high food availability can compensate for the depressing effect of low temperature on modeled larval Arctic cod growth (Thanassekos and Fortier, 2012). Therefore, to get a comprehensive view of the interactive effects of early life traits and environmental conditions on larval Arctic cod growth, including in models food availability and energetic quality in addition to temperature and advection patterns is important.

Here we simulated early life development of Arctic cod under realistic environmental conditions extracted from an ocean and sea-ice model hindcast. For the first time, a bioenergetic model of Arctic cod (David et al., 2022) was coupled to an individual-based Lagrangian tracking model to simulate larval and juvenile growth, specifically accounting for their metabolic processes during dispersal. We tested the sensitivity of modeled prewinter length to variations in early life history traits, such as spawning timing and size-at-hatch, and environmental conditions along dispersal pathways, such as sea-ice cover, temperature, and food quantity. We quantified the recruitment success into juvenile stage before winter (i.e., fraction of larvae reaching a minimum prewinter length) under normal (baseline), low-food, and high-temperature scenarios in four Arctic regions. The present work aims to answer the following questions. Which early life history traits and environmental conditions affect most age-0 Arctic cod prewinter length? How sensitive is larval growth to reduced food availability over the growth season? Will ongoing warming enhance the recruitment success into juvenile stage before winter?

We used an individual-based model (IBM) to simulate temperature- and food-dependent development of eggs, larvae, and juveniles of Arctic cod. The IBM couples a bioenergetic model with a Lagrangian particle tracking model driven by a fully coupled ocean and sea-ice model used to simulate realistic environmental conditions in the Arctic Ocean in which individual eggs, larvae, or juveniles develop and grow while dispersing with ocean currents.

Ocean circulation, water temperature, and sea-ice concentration were generated with the fully validated Biology/Ice/Ocean Modeling and Assimilation System (BIOMAS hereafter; Zhang et al., 2010; 2015). The ocean circulation model in BIOMAS is based on the Parallel Ocean Program (POP) developed at Los Alamos National Laboratory (Smith et al., 1992), nested onto a global ice-ocean coupled model (Zhang, 2005). The ocean model is coupled with a multicategory thickness and enthalpy distribution sea-ice model (Zhang and Hibler, 1997; Zhang and Rothrock, 2003). The coupled ocean and sea-ice model has a generalized orthogonal curvilinear mesh that covers the Northern Hemisphere from 39°N to the geographical North Pole, with the “north pole” of the mesh displaced over the land of Alaska (Zhang et al., 2015). North of 65°N, the horizontal resolution ranges from approximately 4 km in the Bering Strait to 100 km in the Barents Sea (Figure S1). On the vertical dimension, the model has a configuration with 30 layers of varying thicknesses. The ocean mixed layer and euphotic zone are configured with a total of 13 layers in the upper 100 m, of which six vertical layers of 5 m each cover the upper 30 m. To compute the physical advection and early life development of individual fish in the coupled models, ocean temperature and meridional (northward) and zonal (eastward) current velocities were daily- and depth-averaged from the ocean surface to 60 m or to the bottom, whichever was shallower. The BIOMAS model output has been used previously for mesozooplankton tracking and life history development in the Arctic Ocean (e.g., Feng et al., 2016; 2018; Kvile et al., 2018).

The individual-based model of Arctic cod (ArcodIBM) was adapted from an Arctic copepod model (Ji et al., 2012; Feng et al., 2016) by modifying the copepod biological module with the early life development of Arctic cod (David et al., 2022). In ArcodIBM, every individual is represented by a state vector, comprising information on state variables: location (latitude and longitude), life stage (egg, yolk-sac larva, post-yolk feeding larva, juvenile), age (days), length (mm), and weight (g) at each time step. Physical advection of individual fish/particle was calculated every 10 min (i.e., physical time step) and biological variables were updated hourly (i.e., biological time step).

The Lagrangian particle tracking module solves the physical advection equation through a fourth-order Runge-Kutta (RK-4) scheme using daily- and depth-averaged velocities in the upper 60 m from the existing BIOMAS output. All model runs were thus configured in a 2-D horizontal and offline mode. The upper 60-m water column coincides largely with the mixed layer depth in the Arctic Ocean (Peralta-Ferriz and Woodgate, 2015; Hop et al., 2021). The decision to restrict larvae to the upper 60-m layer and model 2-D dispersal was based on the knowledge that larval Arctic cod remain under the sea ice or in the surface open water until late summer (Geoffroy et al., 2016; Levine et al., 2021). The surface distribution of larval Arctic cod was confirmed by a study in the Chukchi Sea by comparing observed larval distributions with modeled larval dispersal under different vertical behaviors (Vestfals et al., 2021). Larval swimming capacity is limited until all fins develop and metamorphosis to juvenile is complete (Ponomarenko, 2000), hence no swimming behavior was simulated in the modeled particles. Given the scale of our study and variability in diffusion terms across the study area, no random-walk was implemented, the focus remaining on the strength of advection in defining dispersal patterns (Largier, 2003). Our choice was based on the knowledge that large-scale oceanographic processes constrain small-scale processes and thus affect larval transport pathways (Pineda, 2000). At model boundaries, no explicit particle behavior was implemented, except a non-absorbing land-ocean boundary condition implemented to allow particles to be advected freely from the land grids by offshore currents. Similar simplifications in representing virtual larvae in IBMs were also considered in other studies when the focus was on large-scale dispersal patterns (Andrello et al., 2013; Álvarez-Noriega et al., 2020) or on rare dispersal events (Treml et al., 2008).

The biological module calculates individual development across four life stages: (1) egg, (2) yolk-sac larva, (3) post-yolk feeding larva, and (4) juvenile (after metamorphosis). Temperature-dependent embryonic development time D (days) was simulated with a laboratory-derived quadratic equation, based on experiments conducted over the temperature range of –0.4°C to 3.8°C (Laurel et al., 2018):

The embryonic development rate ER_t then becomes a reciprocal of development time D, estimated at each biological time step t as a fraction of the development time at a given temperature T.

The embryonic development rate cumulates with each time step until $∑t=1n(ERt)≥1$⁠, when the larva is hatching. The embryonic development duration (days) was estimated as the number of days from spawning until hatching. After hatching, progression from one stage to the next is defined by larval standard length (hereafter referred to as length); hence, a larva with a length <8.5 mm is considered a yolk-sac larva, a larva with length >8.5 mm becomes a post-yolk feeding larva (Michaud et al., 1996), and any fish with standard length >30 mm (reported full range 25–35 mm) becomes a juvenile (Ponomarenko, 2000; Walkusz et al., 2011; Laurel et al., 2017). Development after hatching used the bioenergetic model of Arctic cod described in David et al. (2022), which estimates both weight (g) and length (mm) at each biological time step (hourly). The conversion of weight (W) to standard length (SL) in the model is based on the regression W_(g) = 0.0055 SL_(cm)^3.19, provided by Geoffroy et al. (2016). Temperature- and food-dependent growth rate (GR, g g^–1 day^–1) for larvae and juveniles was computed as the energy gained through food consumption, minus the energy lost through respiration and other metabolic processes, through the general energy-balance equation:

where P is the fraction of the maximum daily consumption C_max, A (= 0.8) is the assimilation efficiency (Hop et al., 1997), representing the remaining consumption after egestion and excretion, R is the basal respiration rate, f_act (= 1.5) is the increase in basal respiration due to activity based on foraging 12 out of 24 h (Lough et al., 2005), and SDA (= 0.375) is the specific dynamic action representing the metabolic increase due to digestion (Hop and Graham, 1995). The equations for consumption (C_max) and respiration (R) are expressed as a function of fish weight and temperature (Text S1; David et al., 2022). In the presence of yolk sac (stage 2), which constitutes the endogenous resources post-hatching, larvae are assumed to consume the maximum daily ration allowed at weight (C_max), therefore the parameter P = 1. After consuming their entire yolk reserves, here using the threshold larval length of >8.5 mm (Michaud et al., 1996), consumption rate becomes a fraction P of the maximum daily consumption for a larva or juvenile (for stages 3 and 4) considering a given prey availability (PA, g_prey g_fish^–1 day^–1), where P = PA/C_max (Roy et al., 2004). When PA < C_max, consumption C becomes a fraction of C_max. The parameter P was estimated previously by fitting the bioenergetic model to multi-year field observations of length-at-age Arctic cod in the Canadian Arctic and the Laptev Sea (David et al., 2022). The value of P per region that best reproduced multi-year length-at-age data was lower in the Laptev Sea (P = 0.7) and higher in Baffin Bay (P = 0.85). In all model runs, a conservative value of P = 0.7 was used for comparability among regions, thus simulating the same feeding conditions in all areas. All modeled and field larval and juvenile length data in this article refer to standard length. Particles found more than 7 days at temperatures higher than 3°C for eggs (Laurel et al., 2018) or higher than 9°C for larvae and juveniles (Geoffroy et al., 2023) were considered dead and eliminated.

Numerical experiments of early life development under realistic environmental conditions were performed in four known spawning areas of Arctic cod, namely the western coastal Chukchi Sea (hereafter Chukchi Sea), northern Baffin Bay (Baffin Bay), waters surrounding Svalbard (Svalbard), and the western coastal Laptev Sea (Laptev Sea; Figure 1). To simulate different monthly spawning events, particles were released on the first day of each month over the spawning period from December to April. While in reality, the release of eggs in the water can take place over the entire spawning period, here we chose for simplification to simulate several spawning events covering the seasonal variability, as commonly seen in other larval fish dispersal studies (Holstein et al., 2014; Petrik et al., 2016; Vestfals et al., 2021). The simulations were conducted for 10 years (2005–2014) to cover the interannual variability. Data for length-at-age of Arctic cod used to parametrize and validate the bioenergetic model originate from the same period (David et al., 2022). Each simulation run ended on September 1 to reflect the fact that Arctic cod ontogenic migration starts at the end of the summer (Geoffroy et al., 2016). An additional comparison run tested September 15 as the ending date. Particles were uniformly released every 0.05° great-circle distance (approximately 5 km) inside the defined spawning areas (Figure 1), resulting in a total number of 5,424 particles released in one spawning event across all four regions. Thus, the number of particles released across all simulated spawning events in the baseline runs summed as 5,424 × 5 months × 10 years = 271,200 in total. Size-at-hatch was fixed at 6 mm, corresponding to the upper range of size-at-hatch observed at low temperatures (Laurel et al., 2018). The high-temperature (High-Temp) scenario was configured for testing the sensitivity of larval growth to an increase in water temperature of 2°C over the growth period (Table 1). The reasoning behind this High-Temp scenario was based on the Paris Agreement (IPCC, 2023) that risks should be minimized if warming is held below 2°C, which should largely support unchanged habitat suitability for Arctic cod. While a monotonous increase in temperature of 2°C over the growth period does not reflect a realistic future warming scenario, our intention was simply to test the response of larval growth to higher temperatures. The reduced food availability (Low-Food) scenario was configured for testing the sensitivity of larval growth to a reduction in food quantity with a decrease of daily ration by 25%, which was mediated through the parameter P in the bioenergetic model (Table 1). The reasoning behind the Low-Food scenario was the uncertainty in the future availability of Calanus copepods, the main prey of age-0 Arctic cod (Bouchard and Fortier, 2020). Hypothetically, even the displacement of Calanus glacialis by the southern, smaller and less lipid-rich species Calanus finmarchicus could lead to decline in Arctic cod (Bouchard and Fortier, 2020), unable to compensate the reduction in the energetic content (Falk-Petersen et al., 2009). The energetic content of prey in the bioenergetic model uses Calanus copepod values (Thanassekos and Fortier, 2012; David et al., 2022). In addition, the size-at-hatch (Hatch-Size) scenario tested the sensitivity of larval growth to variability in size-at-hatch using additional simulation runs in which larval length was initialized with 4, 5, and 7 mm to cover the length-at-hatch variability observed in the field (Sakurai et al., 1998; Ponomarenko, 2000).

For every simulated spawning event with a specific spawning area and a starting year/month, we saved individual fish lengths (mm) on September 1 (hereafter referred to as prewinter length, PWL). To test for significant difference in PWL among scenarios, reduced datasets were used in which PWL was averaged per spawning area, year, month, and scenario. The significant difference in PWL among scenarios was tested with one-way ANOVA, and the homogeneity of variances was tested a priori with Bartlett test. To test interactive effects between scenarios and the factors spawning area, year, and month on PWL, two-way ANOVA was used for testing individual interactions of scenario with each factor.

For each individual fish, the recruitment success into the juvenile stage was then estimated as equal to zero if PWL was smaller than the minimum juvenile threshold length (30 mm) and equal to 1 if PWL was equal to or larger than the threshold. The overall recruitment success in each area per spawning year and month was estimated as the fraction of larvae reaching the threshold prewinter length. The recruitment success was further used as a comparison metric among different scenarios as listed in Table 1. The significant differences in recruitment among scenarios were tested with Kruskal-Wallis non-parametric test. Additionally, we tested the effect of simulation ending time (i.e., the beginning of ontogenic migration on September 1 versus September 15) on PWL and the effect of minimum juvenile length at metamorphosis (25 mm, 30 mm, and 35 mm) on recruitment success.

Environmental conditions experienced by virtual larvae along the advective pathways were recorded daily. For each individual particle or virtual larva, we summarized over their transport pathways the mean water temperature (°C), number of days with sea-ice cover >15% (i.e., the recognized threshold between ice cover and open water), mean bathymetry (m), and dispersal distance (km). Firstly, linear models with Gaussian distributions were used to test the effects of environmental parameters and embryonic development duration on the modeled larval length on September 1. Because temperature was included as a parameter in the bioenergetic model to simulate growth and was shown to explain a large part of modeled length variability (David et al., 2022), we did not use it as a predictor variable of the modeled PWL. We started by testing the significance with each single explanatory variable on PWL, then combinations of two, three, and all four variables. Additionally, we tested the significance of two-factor linear models with interaction terms. For each number of variables, we kept the best model with the highest proportion of explained variance. Secondly, in all the linear model combinations (Table S1), the spawning area, year (2005–2014), and month (December–April) were included as random explanatory variables, tested independently or nested using linear mixed effect (lme) models with the nlme R package (Pinheiro et al., 2020). Lastly, given the limitation of lme to estimate the proportion of the explained variance, the factors with the strongest random effects in the lme models were added as additional predictors to the best linear models previously selected. The model fit and the significance of explanatory variables were evaluated based on the Akaike Information Criterion and the R-squared values. Residual plots were visually inspected for random patterns to assess the adequacy of the models. All the distribution data were checked visually with histograms for normality to satisfy model assumptions, and variables were transformed individually as needed to reach normality (square root: transport distance; cube root: bathymetry). All statistical analyses were performed in R, version 3.6.1 (R Core Team, 2020). Maps were produced using the R package “ggOceanMaps” (Vihtakari, 2024).

Modeled PWL on September 1 ranged from 12 mm to 67 mm over all simulated spawning events, with the regional average ranging from 30 mm in Chukchi Sea to 42 mm in Svalbard (Table 2 and Figure 2). Adding temperature-related mortality did not significantly change the PWL (ANOVA: F_{1, 398} = 0.549, p = 0.459). Intra-annual variability was present in each spawning region without a similar trend across regions (Figure S2). PWL agreed with field observations of length-at-age, estimated from otolith measurements in Baffin Bay and Laptev Sea (Figure S3). Interestingly, the highest PWLs were modeled for Baffin Bay and Svalbard, which had the lowest and highest regionally averaged temperatures over the larval dispersal pathways, respectively. In all simulated spawning months and years, early life stages experienced a regionally averaged sea-ice cover lasting more than 5 months (Table 2). The transport distance was in the range of hundreds of kilometers in all regions, with the distinction that larvae in the Chukchi and Laptev seas were dispersing over the shelf, while in other regions offshore transport prevailed (Table 2 and Figure S4).

Longer ice-covered periods were associated with larger individuals by September, as shown by a significantly positive effect of the number of days with ice cover on PWL (linear model: R² = 0.43, F_{1, 235516} = 179100, p < 0.001; Figure 3 and Table 3). The best two-factor model included the number of days with ice cover in interaction with the embryonic development duration (linear model: R² = 0.66, F_{3, 235514} = 154800, p < 0.001). The duration of embryonic development had a significantly negative effect on PWL, suggesting a shorter development time resulted in larger individuals by the end of the summer (Figure 4). Adding the interaction of the spawning month with embryonic duration significantly improved the model performance indicating a strong effect of the spawning month (R² = 0.89, F_{10, 235507} = 195100, p < 0.001). The effect of the spawning month was also found significant in all linear mixed models, with the strongest random effect of the month nested in area (Tables 3 and S1). The most parsimonious model explaining the largest variability in modeled PWL was represented by the model that combined the interaction terms between the number of days with ice cover and spawning area and the embryonic development duration and spawning month (R² = 0.91, F_{16, 235501} = 147800, p < 0.001).

In the High-Temp scenario, PWL showed an average increase from the baseline scenario of 2.6 mm in Baffin Bay, 2.2 mm in Laptev Sea, and almost 6 mm in Chukchi Sea and a decrease in PWL by 4 mm in Svalbard (Figure 5). Temperature-related mortality increased from 13% in the baseline simulations, all four regions included, to 21% in the High-Temp scenario, with only minor decrease in the average difference in PWL among scenarios (without mortality: PWL increase in Baffin Bay and Laptev Sea by 2.5 mm, Chukchi Sea by 5 mm, and decrease in Svalbard by 4 mm). Because the large sample size might enhance the statistical significance, we tested the significant differences in PWL among scenarios in a reduced dataset, in which PWL was averaged by month, year, region, and scenario. In the High-Temp scenario, PWL did not show a significant difference from the baseline (ANOVA: F_{1, 398} = 1.428, p = 0.233). However, there was a significant effect in PWL with the interaction of scenario with spawning month, and with the interaction of scenario with spawning region, but not with year (Table S2). The Low-Food scenario had a significant effect on PWL (ANOVA: F_{1, 398} = 161.8, p < 0.0001). In the Hatch-Size scenario, size-at-hatch alone had a moderate effect on modeled PWL (ANOVA: F_{3, 476} = 3.806, p = 0.0102; Figure 6).

Spawning events in December, January, and February resulted in annual recruitment rates, defined as the proportion of larvae reaching the juvenile stage, exceeding 0.5 in all regions except Chukchi Sea (Baseline in Figure 7). April spawning events resulted in annual recruitment <0.5 in all regions, with no individual recruited for 3 years in Baffin Bay, 4 years in the Laptev Sea, and 9 years in Chukchi Sea. Recruitment was affected significantly by a decrease in daily food rations in the Low-Food scenario (Kruskal-Wallis: chi-squared = 94.59, df = 1, p < 0.0001). The regions most affected by the decrease in food rations were Chukchi and Laptev seas where recruitment declined 5-fold and 3-fold, respectively (Low-Food scenario in Table 4). March and April spawning events were affected strongly in all regions by the decrease in food quantity, resulting in zero annual recruitment for most of March and for all of April events (Figure 7). Alternative hatch sizes (tested sizes: 4 mm, 5 mm, and 7 mm) induced a small, yet a significantly positive effect on the overall recruitment success (Kruskal-Wallis: chi-squared = 15.389, df = 3, p = 0.0015; Tables 4 and S3). In the High-Temp scenario, early spawning events (December–February) maximized recruitment in all regions while the latest spawning event in April minimized recruitment, with nearly all to zero. The spawning event in March maintained comparable recruitment rates with the Baseline scenario across regions (Figure 7). In spite of opposite monthly trends in recruitment with the High-Temp scenario, overall there was a significant effect of increased temperature on recruitment (Kruskal-Wallis: chi-squared = 3.5204, df = 3, p = 0.0060).

Recruitment changed linearly with the juvenile threshold length at metamorphosis (tested lengths of 25 mm, 30 mm, and 35 mm; Table 5) decreasing on average from 0.86 for the 25 mm threshold to 0.46 for the 35 mm threshold. Sampling 2 weeks later would result in an increased recruitment by 0.11–0.16 over all threshold lengths (Table 5).

Arctic cod early life stages are known to develop slowly, often spending 5–6 months at sub-zero temperatures under sea ice before using the ice-free season to boost their growth before winter sets in (Bouchard et al., 2017; LeBlanc et al., 2020). Reaching large prewinter length and storing sufficient amounts of lipids increase winter survival for age-0 Arctic cod (Bouchard and Fortier, 2020; Copeman et al., 2020; 2022). Our modeled PWL confirms that winter spawns result in the largest individuals by the end of summer and that the PWL range results from an extended spawning season rather than environmental variability. The long spawning season and iteroparous (repeated) spawner strategy of Arctic cod (Hop et al., 1995) are related to the adaptive trait of Arctic cod to reproduce during the polar night, giving Arctic cod an advantage over boreal fish species (Geoffroy and Priou, 2020). This adaptive trait ensures that environmental risk of extended periods of low temperatures for early life stage development and mismatch of hatching with spring production is spread over multiple months. The reproductive strategy of Arctic cod is thus believed to maximize the PWL, leading to increased survival of the age-0 group and juvenile recruitment (Bouchard and Fortier, 2011; Bouchard et al., 2017).

Mimicking the long spawning season from December to April over 10 years, our model captured well PWL variability in four Arctic regions, that is, Baffin Bay, Chukchi Sea, Laptev Sea, and Svalbard, indicating significant differences by spawning year. Yearly differences in larval and juvenile lengths were also reported at the end of summer in Chukchi Sea (Vestfals et al., 2019; Deary et al., 2021), Svalbard (Eriksen et al., 2015; 2020), Baffin Bay, and Laptev Sea (Bouchard and Fortier, 2008; 2011). Cold versus warm years explained yearly differences in PWL in Chukchi Sea through temperature variations and modified circulation patterns affecting larval transport, resulting in different growth conditions (Vestfals et al., 2021), albeit the classifications of cold and warm years were not the same across all Arctic regions (Alabia et al., 2023).

In comparison to spawning timing, variations in size-at-hatch induced only small effects on PWL, according to our analysis. Variation in size-at-hatch is often attributed to the effect of hatch timing within a batch of eggs or to population differences (Sakurai et al., 1998; Laurel et al., 2008). In laboratory, delayed hatching under lower temperature resulted in Arctic cod being 1.5 mm larger at hatch at –0.5°C than at 4°C (Laurel et al., 2018). The length-temperature relationship in wild Arctic cod larvae is not straightforward, given the large range of size-at-hatch (4–7 mm) observed at sub-zero temperatures, thus population differences could be presumed to be the cause of this variability (Sakurai et al., 1998). Given the minimal effect of size-at-hatch on PWL, using an average size in our simulations was considered reasonable for computational constraints without losing information.

Probably the most important trait determinant of early life stage development is spawning-ground selection by adults, initiating the dispersal patterns and hence predetermining the environmental conditions encountered by eggs and larvae during dispersal (Figure 8; Eriksen et al., 2020; Gjøsæter et al., 2020). However, the exact location of spawning grounds remains a major unknown in the ecology of Arctic cod, due to lack of direct observations in winter, with a few exceptions mentioned in old Russian literature (Aune et al., 2021) or recently informed by Inuit knowledge (Bouchard et al., 2023). Spawning grounds in Svalbard and Chukchi Sea were inferred from the age-0 Arctic cod distribution at the end of summer matching best with simulated particle dispersal (Vestfals et al., 2019; Eriksen et al., 2020; Vestfals et al., 2021). In Baffin Bay and Laptev Sea, no similar work has been performed, so spawning grounds could only be presumed from summer larval distributions and expert knowledge. Therefore, the particle release area in our simulations should be regarded rather as potential spawning habitats within the two regions. In spite of uncertainty in spawning grounds, modeled PWL reproduced well field PWL variability within and between regions. Larger PWL was reported in Svalbard (40 mm mean) (Eriksen et al., 2015; 2020) and Baffin Bay (Bouchard and Fortier, 2011) compared to Chukchi (30 mm mean; Vestfals et al., 2019; Deary et al., 2021) and Laptev Sea (26 mm and 35 mm, for early versus late hatchers, respectively; Bouchard and Fortier, 2008). Larger PWL in Svalbard could be explained by the inflow of warm Atlantic water (Huserbråten et al., 2019), while the smaller PWL on the eastern Chukchi Sea shelf, even though located at a more southern latitude, could result from the influence of the cold coastal Siberian current (Vestfals et al., 2021). Low temperatures on the Siberian shelf prevail in the Laptev Sea and maintain slow growth in Arctic cod early life stages there (Bouchard and Fortier, 2008; 2011). In spite of generally low temperatures in north Baffin Bay, early warming of surface water following spring polynyas is believed to induce faster growth in this region, resulting in larger PWL of wild Arctic cod (Bouchard and Fortier, 2011), similarly reflected in our modeled PWL. PWL seems to largely reflect regional sea-ice dynamics and temperature variations within these Arctic regions. Nonetheless, while modeled PWL captured well regional and temporal variations of length-at-age in the field, the dominant life history trait in explaining the largest PWL variability is spawning month, suggesting the significance of an extended spawning season.

Life history traits in fish are inherited, but environmental variability acts upon these traits either through expected seasonal cycles or stochastic yearly or regional events (Hutchings, 2021), events that in the Arctic induce changes in sea-ice dynamics and temperature among others (Peralta-Ferriz and Woodgate, 2015). Among the environmental variables recorded along the dispersal pathways of eggs and larvae, we found the duration of sea-ice cover to be a significant parameter for PWL, with more days with ice cover resulting in larger PWL. The positive effect of sea-ice cover could be related to spawning month, as winter spawns experience more days with ice cover than spring spawns.

In addition to the number of days with sea-ice cover, the largest part of modeled PWL variability was explained by the embryonic development duration and the spawning month. Our findings suggest that larvae from early spawning events which experience shorter embryonic duration, hence hatch earlier, result in largest PWL, likely due to an extended growth season. A positive correlation between hatching timing and ice breakup previously has been shown to determine larger length-at-age in larval Arctic cod, highlighting a positive adaptation of this early life history trait to environmental variability (Ponomarenko, 2000; Bouchard and Fortier, 2011; Aune et al., 2021; Herbig et al., 2023). In Baffin Bay, the hatching peak period was estimated in April–May (Bouchard and Fortier, 2011). Compared to the Canadian Arctic, hatching in Svalbard was observed starting as early as March and extending to late June (Eriksen et al., 2015), while in Laptev Sea hatching can extend from January to July (Bouchard and Fortier, 2008). A long growth season post-hatching overcompensates for the slow growth induced by low temperatures during the first months of embryonic development and ultimately results in larger PWL (Bouchard and Fortier, 2011), but decreases survival of early hatchers due to longer exposure to predation (Fortier et al., 2006; Gjøsæter et al., 2020). We did not consider mortality due to predation in our simulations, as our focus was on the effect of environmental variability on larval growth. Moreover, because we reported monthly averages in modeled PWL, our results remain unaffected by the non-inclusion of predation. We did include a temperature-related mortality of eggs and larvae, considering upper thermal limits of 3°C for eggs (Laurel et al., 2018) and 9°C for larvae (Geoffroy et al., 2023), which eliminated mostly outliers in our data (Figure S5) and, consequently, did not significantly affect the modeled PWL.

Our choice of study regions was based on location variability (inflow of Atlantic or Pacific water, Canadian archipelago, and Siberian shelf) and availability of field data for validation and comparison. Interestingly, across these four regions, the life history traits of Arctic cod, specifically the embryonic duration and the spawning timing, rather than the region or environmental variability, such as the number of days with ice cover, had the largest effect on PWL. These life traits seem to dictate primarily the success of a large fraction of larvae reaching the juvenile stage before winter. Under current environmental conditions, the selection of such specific traits might be favoring the dominance of Arctic cod in the Arctic Ocean.

Warming over the past three decades suggests a positive effect on age-0 Arctic cod growth, indicating an increase in length-at-age (Dupont et al., 2020). Our High-Temp scenario also suggests a positive effect of increased temperature on PWL, at least in three of the four regions studied. While optimal temperatures for larval Arctic cod growth are reported to fall between 2°C and 5°C (Koenker et al., 2018b; Eriksen et al., 2020), larvae are known to dwell under the sea ice at sub-zero temperatures during a large part of their development. Our simulations also showed that larvae spend several months under the sea ice. Being well adapted at near-zero temperatures, larval Arctic cod have for now competitive growth and survival advantage over sub-Arctic fish species (Koenker et al., 2018a; 2018b; Laurel et al., 2018; Vestfals et al., 2019). Moreover, sub-Arctic fish larvae do not seem to compete with Arctic cod for the same food resources (Falardeau et al., 2014; Bouchard et al., 2022), but with ongoing warming, the balance might shift in favor of sub-Arctic species (Renaud et al., 2012; Christiansen, 2017; Koenker et al., 2018a; Geoffroy et al., 2023). Under a low-food scenario tested in the laboratory, Arctic cod maintained competitive growth at low temperatures, but they lost the advantage at higher temperatures (Koenker et al., 2018b). Larval and juvenile Arctic cod are physiologically capable of compensating for an increase in temperature with higher food intake, if enough prey is available (Kunz et al., 2016). In their natural habitat, however, acquiring more prey implies an increased foraging effort, making the net energy gain uncertain. To test the consequences of a low-food scenario in the natural habitat of Arctic cod, we simulated a decrease in larval daily ration by 25%, a value implying slower, yet positive growth, over the range of temperatures encountered during simulated dispersal pathways (David et al., 2022). A decrease of 25% in daily ration would halve the recruitment of the age-0 group in September, which could have dire implications for the juvenile stock. Even with less food available, larvae originating from early spawning events (December–February) still had high chances of reaching metamorphosis length by the end of summer. The risk of mismatch with zooplankton production is indeed higher for early hatchers, but slow metabolism at sub-zero temperatures and delayed yolk-sac exhaustion would allow Arctic cod larvae to sustain food deprivation for weeks until food becomes available (Copeman et al., 2022; Koenker et al., 2018a), a decisive advantage over sub-Arctic species. For late hatchers, our results show that food deprivation could have immediate consequences, leading to recruitment failure.

With high prey availability and slightly higher water temperature, larvae could reach the juvenile stage rapidly, increasing swimming and feeding capabilities, thus enhancing their survival (LeBlanc et al., 2020). However, should the observed ice reduction and warming continue, juvenile recruitment collapse is anticipated (Huserbråten et al., 2019; Bouchard et al., 2021). According to the Paris Agreement (IPCC, 2023), risks should be minimized if warming is kept below 2°C, which should largely support unchanged habitat suitability for Arctic cod. In Baffin Bay (≥55°N), from 2005 to 2100, median ocean temperatures are projected by IPCC scenario RPC8.5 to increase by 1.9°C in the surface layer during spawning season December–March (Cote et al., 2021). As shown by our results, Arctic cod in northern Baffin Bay would benefit from a 2°C temperature increase by speeding up embryonic development and larval growth, resulting in larger PWL. However, the temperature-related mortality increased from 13% to 21% over the growth period. The suitable habitat for Arctic cod egg survival, on the other hand, at least in the southern parts of Baffin Bay, will be restricted spatially in this case and shifted northward (Cote et al., 2021). In warmer areas around Svalbard, such a temperature rise will have a detrimental effect on larval growth, as suggested by our results indicating a decrease in average PWL by 4 mm. The decrease in PWL is due to the dome-shaped growth rate with temperature defined in the bioenergetic model, where growth rate decreases above the thermal optimum, limiting the increase in length (David et al., 2022). A limitation of the bioenergetic model resides in using data from the high Arctic to parametrize the temperature-dependent function of food consumption (Thanassekos and Fortier, 2012), in the absence of data covering the entire temperature range of larval Arctic cod habitat. With an updated function using extended temperature data, modeled growth rate might capture better Arctic cod response and possible adaptation to warming.

Ongoing warming, delays in winter onset, and hence longer growth season are expected in the Arctic Ocean (AMAP, 2021). Our results show that extending the growth season by 2 weeks would allow Arctic cod to reach a larger PWL before winter onset and increase the fraction of larvae metamorphosizing into juveniles. Warming could also potentially induce early metamorphosis to the juvenile stage, at a shorter length (Benoît et al., 2000); thus, considering a shorter threshold PWL of 25 mm for passing metamorphosis compared to 35 mm (i.e., the extremes of the metamorphosis length range in Arctic cod) would double the recruitment to the juvenile stage in September (Figure 8). Therefore, if sufficient food is available, current conditions of early ice retreat, slight warming, and delay in ice freeze-up should reflect an increase in the fraction of larvae metamorphosizing into juvenile stage, hence in recruitment success, in September. Nonetheless, changes in size at metamorphosis influence predator avoidance, prey capture, and competitive abilities. Therefore, not only changes in the length-at-age relationship with warming should be of interest in understanding the development of early life stages of Arctic cod, but also changes in length- and age-at-metamorphosis because such changes would also affect recruitment.

We demonstrated that early life history traits significantly affect larval length before winter and could therefore determine recruitment success to the juvenile stage. Environmental conditions such as temperature and food availability indeed mediate larval growth, but ultimately a long growing season, provided by early spawns and short embryonic duration, seems most beneficial for larvae to reach the metamorphosis length to juvenile before winter. Nonetheless, the recruitment success to juvenile stage comprises the overall value of different interacting factors coding life history traits and environmental forcing encountered by early life stages during dispersal. We hypothesize that our estimation of recruitment success using a threshold length at metamorphosis represents an effective indicator of the age-0 group recruitment to the pelagic juvenile population. Validating this metric with extensive multi-year data would be a crucial next step toward predicting age-0 recruitment success, thus to significantly advance our prediction of future changes in Arctic cod populations.

The R scripts generated in this study for data analysis and figure plotting are publicly available on Github and can be accessed at https://github.com/cld321/ArcodIBM. The model outputs and the scripts for post-processing are stored on the WHOI server at https://ulysse3.whoi.edu/thredds/catalog/data/ArcCodIBM_output/catalog.html.

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

Supplemental Material (PDF)

This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund, through the Ocean Frontier Institute who supported CLD. ZF was supported by the National Natural Science Foundation of China (Grant# 42176225). IDA was supported by Japan’s national program on Arctic studies, Arctic Challenge for Sustainability II (Grant number: JPMXD1420318865).

The authors declare no conflicts of interest.

Contributed to conception and design: CLD, RJ, and JAH.

Contributed to acquisition of data: CLD, CB, ZF, and JZ.

Contributed to analysis and interpretation of data: CLD, RJ, JAH, CB, and HH.

Drafted and/or revised the article: CLD, RJ, JAH, CB, HH, IDA, ZF, and JZ.

Approved the submitted version for publication: CLD, RJ, CB, HH, IDA, ZF, and JZ.

How to cite this article: David, CL, Hutchings, JA, Feng, Z, Bouchard, C, Alabia, ID, Hop, H, Zhang, J, Ji, R. 2025. Effects of early life history traits and warming on Arctic cod prewinter length and recruitment. Elementa: Science of the Anthropocene 13(1). DOI: https://doi.org/10.1525/elementa.2024.00015

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

Associate Editor: Maxime Geoffroy, Fisheries and Marine Institute of Memorial University of Newfoundland, St. John’s, Newfoundland, Canada

Knowledge Domain: Ocean Science