The interactive effects of temperature and food consumption on growth of larval Arctic cod (Boreogadus saida): A bioenergetic model

Understanding the factors determining Arctic cod (Boreogadus saida) recruitment is a key element in predicting the response of the species to climate change in different parts of its circumpolar distribution. The Arctic cod is integral to Arctic ecosystems because of its high biomass and its importance as a food source for seabirds, fish, and marine mammals, especially in the High Arctic (Mecklenburg and Steinke, 2015; Mueter et al., 2020). The life cycle of Arctic cod is intimately connected to sea ice, suggesting a particular vulnerability of this species to changes in sea ice, especially during its early life stages (Graham and Hop, 1995; Bouchard et al., 2017). Larval growth rate could directly influence survival probability during early life and, thus, recruitment to the adult population. Temperature and food consumption are the most important extrinsic factors of Arctic cod larval growth (Thanassekos et al., 2012; Koenker et al., 2018), but determining their relative importance and their interactive effects remains a great challenge, especially for wild larvae.

The reproductive strategy of Arctic cod includes a very long hatching season, extending from January to July (Bouchard and Fortier, 2008, 2011). Larvae from the same region and year but with different hatch dates can experience very different environmental conditions during their development, resulting in different growth and survival rates (Bouchard and Fortier, 2008, 2011). In region–years with environmental conditions allowing for the survival of winter hatchers, favored by increased food availability and lower predation, these individuals have an important advantage over the summer hatchers (Bouchard et al., 2017; LeBlanc et al., 2020). Despite slower growth during the first weeks of life resulting from low temperature, a longer growth season maximizes their prewinter weight and lipid reserves, important assets to surviving the first winter (Sogard, 1997; Copeman et al., 2020).

Temperature and food availability in the natural habitat affect larval Arctic cod growth by regulating their bioenergetic balance between the food intake and metabolic losses. However, the interaction or combined effects of temperature and food on growth remain unclear. Studies on larval fish at low temperatures have found that growth rates documented in the laboratory can be higher than those observed in the wild (see examples in Clarke and North, 1991; Icelandic Atlantic cod Gadus morhua in Björnsson, 1999). A key difference between laboratory and natural settings could be the food environment in terms of both quantity and quality. If optimal food conditions are not met by larval Arctic cod in the wild, then answering the question of whether temperature imposes an upper limit on the growth from the spatial distribution data alone is difficult (Clarke and North, 1991; Eriksen et al., 2020). Laboratory measurements of metabolic rates, such as respiration, heart rate, and feeding, suggest that larval Arctic cod exhibits a broader thermal tolerance than what is realized in their natural habitat (Drost et al., 2016; Kent et al., 2016). Therefore, to consider temperature as the determining factor of their natural habitat range is too simplistic. The challenge resides in considering the limitations imposed by both temperature and food on metabolic rates in order to be able to model larval growth at comparable levels to what is observed in the wild.

For larval Arctic cod, the relative sensitivity of growth to temperature and food consumption has not been directly compared over the full range of values observed in the field. A modeling study by Thanassekos and Fortier (2012), which used a bioenergetic model to simulate larval growth at different temperatures and food levels, predicted growth rates within the observed range in the field, but only for larvae smaller than 12 mm. Recent experimental work on early larval development in relation to temperature has revealed a growth optimum of 3°C–6°C (Koenker et al., 2018), surprisingly higher than estimated optimal temperatures derived from field observation data (Bouchard and Fortier, 2011). While the higher temperature optimum for growth might hold true for larger larvae (i.e., ≥12 mm) and juveniles, previous studies reported that temperatures above 3.5ºC result in poor hatching success and larval deformities (Sakurai et al., 1998; Kent et al., 2016). Moreover, direct comparisons between any kind of model predictions of growth patterns and observed growth of larvae in natural or seminatural environments are lacking for larval Arctic cod.

In the present study, we adapt and expand the bioenergetic model originally developed by Thanassekos and Fortier (2012) to cover the full-size range of larval Arctic cod and to test its applicability over different High Arctic regions. The model is used as an integrative tool for investigating the effects of temperature and food consumption on Arctic cod larval growth in order to understand the differences between experimental work and field observations. We hypothesize that (1) variations of 1°C–2°C at the lower temperature range can lead to large differences in size-at-age of larval Arctic cod by the end of their first summer and that (2) both food quantity, in terms of fraction of maximum daily ration consumed, and food quality, in terms of energetic content of prey, are equally important in driving the observed variability of growth rates. This study will help us understand how the Arctic cod larval growth is affected by the interactions of early life history traits (e.g., hatching timing and size-at-hatch) and environmental conditions (e.g., temperature and food availability), which represents an important step in predicting how this Arctic species will respond to ongoing environmental changes.

Previously published length-at-age data of age-0 Arctic cod were estimated from otolith measurements of pelagic larvae and juveniles sampled in late summer and early fall between 2003 and 2015 in four regions of the Arctic Ocean (Bouchard and Fortier, 2011, 2020), for a total of 14 region–years. Fish were sampled from August 14 to October 8, standard length (SL) of fresh individuals measured, and otoliths extracted. A total of 1,092 fish were aged by otolith analysis (Table 1).

The sampling was conducted in three regions of the North American Arctic (Baffin Bay, Beaufort Sea, and Kitikmeot region) and one region of the Siberian Arctic (Laptev and East Siberian Sea, hereafter referred to as Laptev Sea). All sampled regions are located north of the Arctic Circle and ice-covered most of the year, with some spatiotemporal differences in timing of sea-ice retreat, ranging from May to September. Daily sea surface temperature (SST) in these regions ranges from near freezing −1.8°C immediately under ice cover to about a maximum 5°C in summer (Figures 1 and S1–S4).

Daily SST fields were downloaded from the Estimating the Circulation and Climate of the Ocean Consortium (https://www.ecco-group.org/home.cgi; ECCO2 cube92 SST data set; Menemenlis et al., 2008) in each region–year, covering the estimated months of hatching and larval growth (from February 1 until September 15). The four regions were delimited geographically to include all Arctic cod sampling stations, as described in Bouchard and Fortier (2011) and LeBlanc et al. (2020).

Bioenergetic models compute growth rate as a balance between food consumption (energy gain) and metabolic losses (respiration, excretion, and other losses through activity), both of which are affected by fish size and temperature (Brett, 1979; Railsback and Rose, 1999). These models can be used to synthesize metabolic processes and infer ecological relationships between fish and their environment (Ricker, 1979; Kitchell and Breck, 1980; Roy et al., 2004).

Here, we adapted a bioenergetic model (BEM hereafter) from Thanassekos and Fortier (2012; TF2012 model hereafter) to simulate Arctic cod larval growth over time and to quantify the effect of temperature and food consumption on larval growth. The model assumes that (1) daily consumption is a constant proportion of maximum consumption, (2) egestion and excretion are constant fractions of consumption, (3) energy content of fish and prey are constant, and (4) the ratio between active metabolic rate and basal standard metabolic rate is fixed. A major difference between our model and the TF2012 model is in the estimation of food consumption. First, a simplified formula, rather than a full forage module, was used in our estimation of consumption due to the lack of experimental data and associated parameter uncertainties. The simplified formula is similar to one used by Kitchell and Breck (1980) and Roy et al. (2004). The simplified consumption estimation in our model does not require input from field observations of prey data, which are not available for the entire growth season in large parts of the Arctic. Second, we revised the fraction of nauplii in the diet used in the TF2012 model, which was estimated for larvae <12 mm in length, and extended the formulation to the full age-0 length range (4–50 mm). To ensure the applicability of our model at the pan-Arctic scale, we validated our results with length-at-age estimations from otolith measurements from four High Arctic regions covering multiple sampling years (Table 1). The model predictions remain restricted to the High Arctic until further validation can be made with data on larval Arctic cod from more southern and warmer geographic areas.

Throughout this article, the length of fish refers to the fresh standard length (SL, mm) unless specifically stated otherwise. The conversion of wet weight (W) to SL in the model is based on the regression W_(g) = 0.0055 SL_(cm)^3.19, provided by Geoffroy et al. (2016).

The weight-specific growth rate (G, g g^–1 d^–1), here defined as daily weight change per unit of fish weight, can be calculated based on the following energy-budget equation:

where C (g g^–1 d^–1) is the weight-specific consumption rate and A is the assimilation efficiency, representing the proportion of consumption remaining after egestion and excretion (Hop et al., 1997). R (g g^–1 d^–1) is the weight-specific rate of respiration and basal standard metabolism, f_act is a proportional increase in respiration due to activity, and SDA is the specific dynamic action. These terms are explained in detail in the sections below.

The weight-specific consumption rate (C, g g^–1 d^–1) is given by the formula:

where C_max is the maximum weight-specific consumption rate at weight (W) and P is defined as a proportion (ranging from 0 to 1) of C_max. The maximum weight-specific consumption rate at weight (C_max, g g^–1 d^–1) is defined as:

where a_c = 0.0337 and b_c = –0.2862 were estimated from the maximum daily ration formula:

experimentally determined for a size range of fish from juvenile to adult (weighing 3.99–59.19 g) at 0°C by Hop et al. (1997). The initial formula log₁₀ MDR (mg_prey d^–1) was written in its exponential form, converted to grams and divided by fish weight (W) to obtain the formulation MDR $(gpreygfish−1d−1)=ac⋅Wbc$⁠, here corresponding to consumption as a function of fish weight.

The temperature function for the consumption f(T) was estimated by Thanassekos and Fortier (2012) from gut content volume (V) of field-captured Arctic cod larvae. The functional response was fitted as the regression of V/V_max to temperature and standardized to 1 (where V_max is the maximum prey volume at larval length at 1-mm intervals), resulting in:

The ratio of calorific density of prey (CD_prey) to calorific density of larvae (CD_fish) was added to express the efficiency of energy conversion from prey to fish tissue (Thanassekos and Fortier, 2012). By adding this ratio, the consumption becomes an expression of the proportional daily increase in fish weight (g_fish g_fish^–1 d^–1), eliminating the energy units. The calorific density of prey changes as larger larva selectively feed on larger prey, gradually switching from a mixed diet of copepod nauplii (CD_n) and eggs (CD_e) to copepod nauplii alone:

The fraction of nauplii (npl) in the diet was estimated from gut content of larvae <12 mm by Thanassekos and Fortier (2012; Table 2; Figure 2). In addition to nauplii, larvae >15 mm gradually include copepodites in their diet; however, nauplii remain the dominant items in the diet of larvae <25 mm (Bouchard and Fortier, 2020). As both groups have similar calorific density, we did not specifically include the copepodites but assume that they are part of the diet as larvae increase in size and improve their swimming ability. Gradual increase in swimming speed with the size of both larvae and a larger prey would result in little difference in the outcome of any encounter as long as a size ratio of 1/10 is maintained (Fiksen and MacKenzie, 2002). Hence, the balance between energetic expenses of foraging, including failed attempts on larger prey, and energetic gain on captured prey would presumably be maintained.

The parameter P in Equation 2 was defined by Roy et al. (2004) as the proportion of C_max possible for a larva at a given temperature and prey availability (PA, g_prey g_fish^–1 d^–1), where P = PA/C_max. When PA < C_max, consumption C becomes a fraction of C_max. Increased PA leads to higher consumption until P reaches 1, as further increases cannot be ingested due to physiological constraints. For yolk-sac larvae, P is equal to 1, as yolk reserves constitute endogenous resources for feeding, hence maximum daily consumption is assumed. Because the conversion of caloric content from prey to fish biomass does not have relevance for yolk-sac larvae, which rely on endogenous feeding, the term CD_prey/CD_fish was excluded from the C_max formulation (Equation 3) for the duration of the yolk sac. When Arctic cod larvae reach a length of >8.5 mm, they are considered to have entirely lost their yolk reserves (Michaud et al., 1996) and to rely exclusively on exogenous feeding, with consumption rate computed as a fraction P of the maximum consumption C_max. We use this simplification in the BEM by replacing a complicated foraging module with a single parameter P. Treated as an unknown, in the absence of data on PA, the parameter P can be determined by fitting the model to field observations.

The weight-specific standard respiration rate (R, g g^–1 d^–1) is modeled as a function of weight and temperature, adjusted by the ratio of the oxycalorific coefficient (Q_ox) to calorific density of fish larva (CD_fish):

where parameters a_v = 0.0355 and b_v = –0.1699 were estimated by combining a weight-specific respiration equation with a temperature-dependent respiration equation, the latter becoming the temperature function for the respiration f(T) in Equation 7. For f(T) formulation, oxygen consumption rates, reported as VO₂ (gO₂ g_fish^–1 h^–1), were shown to increase linearly with temperature:

The oxygen consumption rates were measured for 5-month acclimated fish at temperatures of 0.4°C–2.7°C and standardized to a 10-g fish for comparison (Hop and Graham, 1995), using the following formula:

in which the weight exponent B = 0.8301 originates from the weight-specific equation estimated by Holeton (1974) at a constant subzero temperature of –1.5°C for Arctic cod in the weight range of 0.6–122.0 g:

Solving Equation 9, in which VO_{2 (10)} is replaced with f(T) expression as in Equation 8, VO_{2 (W)} (gO₂ g^–1 h^–1) can be written as a function of both weight and temperature and becomes the R expression in Equation 7. To this formulation, the oxycalorific coefficient (Q_ox, cal g O₂^–1) was added to convert weight units of O₂ into calories, which divided with the calorific density of fish larva (CD_fish) will translate respiration into proportional daily loss in weight (g_fish g_fish^–1 d^–1), thus using the same units as consumption and solving the main energy-budget equation (Equation 1). The oxycalorific coefficient was estimated considering the respiration substrate for feeding Arctic cod with a diet mainly on Calanus copepods (Hop and Graham, 1995, and references therein).

Because the weight-specific respiration was estimated for resting fish, an additional increase is required to reproduce an active metabolism. When larvae are foraging, the standard metabolism is elevated, estimated by a factor of 2.0 for larval Atlantic cod (Lough et al., 2005); assuming that Arctic cod larvae are foraging during the daytime (12 of 24 h), the standard respiration rate was increased by a factor of f_act = 1.5. The additional energy expended on all activities of the body incidental to the ingestion, digestion, absorption, and assimilation of a meal was quantified by the SDA and was expressed as daily proportional increase in respiration by Hop and Graham (1995), where SDA = 0.375.

To find the best value of the parameter P, expressing the proportion of the maximum daily food consumption, length-at-age data from otolith measurements from four Arctic regions were used. Each of the individual region–year data sets shown in Table 1 (we used 10 of the 14, as the year 2005 was not used in parameter fitting but kept for a second independent validation) was randomly split in two (60%, 40%) for either fitting parameter P or for model validation, respectively. The choice of 60%–40% was made gradually to increase the statistical power for fitting P of some limited data sets, which had only few data points. Improving the fitting of the parameter P ultimately led to an increase in modeling efficiency (EF), as explained in the next paragraph. For each region–year, the model was initialized with the minimum larval length sampled in that respective region–year and run for the number of days to cover the growth of the oldest larvae (i.e., age in days of the oldest larvae minus the age of the youngest). The model used regionally averaged temperature (as showed in Table 1), and a range of P values from 0.5 to 0.9 (at 0.01 intervals) were tested.

For each value of P tested, an EF metric was calculated based on the formula EF = 1 – Σ|y_o – y_p|²/ Σ|y_o – y_m|² (Mayer and Butler, 1993), where y_o are the length-at-age values observed, y_p are the values predicted, and y_m is the mean of observation values. The P value with the highest EF was selected for the respective region–year. Then, using this value of P, the model was run again, and EF was calculated, using the second subset kept for model validation, representing 40% of the data in that respective region–year. Due to high uncertainty in the daily location of the sampled larvae, which in the High Arctic is on the order of a few hundred kilometers over the growth period (Eriksen et al., 2020; Vestfals et al., 2021), the choice was made to use regionally averaged temperature for model simulations. Moreover, different hatching cohorts can also experience different temperatures, such as early hatcher versus late summer hatcher. To override the growth differences arising from different hatching cohorts, the second choice made in this regard was to use an average temperature over the growth period of the dominant cohort, that is, the April hatchers (mean SST values between April 15 and September 15, shown in Table 1).

We acknowledge the potential errors associated with using an average temperature to simulate growth over few months. Therefore, we used an additional data set, which is the year 2005 in all four High Arctic regions (four region–year data sets), to evaluate the model efficiency in predicting growth for different hatching cohorts using daily temperatures. For this validation, an average value of P was used per region and growth was initialized with four larval lengths, 4, 5, 6, and 7 mm, covering the hatching size variability of larval Arctic cod (Ponomarenko, 2000). The model was run using regionally averaged daily temperatures from the ECCO2 model (Menemenlis et al., 2008) and simulating larval growth with different daily hatching dates, starting from February 15. All simulations ended on September 15. EF was calculated for each region, using the predicted lengths within the same age interval as observations.

The model sensitivity analysis was conducted to assess the relative importance of parameters to the modeled growth rates. Here, we focus the analysis on the temperature and the parameters related to food quantity (P) and quality (i.e., CD_n, CD_e, CD_f). Sensitivity analysis for the full list of parameters can be found in Supplemental Material (Texts S2 and S3). First, we quantified the sensitivity of growth rate to variation in model parameters through individual parameter perturbations (IPP; Bartell et al., 1986). Each parameter value was varied ±20% of the baseline value at a time while the other parameters were held constant (Table 2). For temperature to be representative of the High Arctic regions, a reference temperature of 0°C was chosen which then was varied ±1°C. Variation in larval length after 30 days of growth, initialized with a size at hatch of 6 mm, was compared to a reference larval length simulated using the nominal values of all parameters. The relative change in modeled length (L_m) from the reference (L_ref) was estimated as the ratio ▵L = (L_m – L_ref)/L_ref. Second, we applied a global sensitivity analysis in which all parameters were varied simultaneously, each randomly selected from a normal distribution. We implemented this method by using Latin Hypercube Sampling (van Griensven et al., 2006) and determined the relative contribution of each parameter variation to prediction error (Bartell et al., 1986). Normal distributions were assigned to the parameters included in the IPP analysis by using the nominal values as the mean and their ±20% variation as the minimum and maximum values (Table 2). Values were drawn randomly from the distributions to create a set of parameter values in each model simulation. Each parameter distribution was sampled 300 times, resulting in a data set 300 × 6 (five parameters and modeled larval length). The data set was analyzed statistically by calculating Pearson correlations between simulated larval length and each parameter. The relative contribution of each parameter to prediction error was estimated by multiple linear regression and a bootstrapping method. We repeated the analysis for two sizes of larva (SL of 6 and12 mm).

Different simulation scenarios were cross-tested in which (1) food quantity and (2) food quality were modified. The change in food quantity was simulated by using a fraction of maximum daily food consumption of P = 0.6 which then was reduced to P = 0.3. The alteration in food quality was implemented in the model by using either a diet based on a Calanus nauplii caloric content of 1716 cal g^–1 (Table 2) or a diet based on Artemia nauplii of 1,295 cal g^–1 (equivalent to 20.9 kJ g^–1 according to Vanhaecke et al. (1983), considering 1 cal = 4.1868 J and a water content of 74%). In all scenarios, the BEM was initialized with a length-at-hatch of 6 mm and used a constant temperature of 0°C during growth. The resultant modeled length-at-age was compared among different scenarios, informed by both field measurements and laboratory experiments.

Within the same region, the best values estimated for the parameter P were in good agreement across the years (Table 3). Higher average P per region of 0.85 and 0.82 were estimated for Baffin Bay and Kitikmeot region, the two regions with the coldest temperatures, followed in decreasing order by Laptev Sea with 0.70 and Beaufort Sea with 0.60. The implementation of these values of P in the validation runs resulted in an EF between 0.6 and 0.97, with the lowest values associated with the region–years having the lower number of observations, irrespective of region (Table 3). This association was also noticeable in the visual comparison of observations versus predicted values of length-at-age (Text S1, Figures S5–S14).

Modeled length-at-age estimated using a fixed value of P per region and a constant temperature, averaged over the years, was well within the range of observations based on otolith measurement from four Arctic regions, when compared to individual regressions fitted per region and year (Figure 3).

When daily temperatures were used to simulate growth of larvae with different hatching dates, starting from February 1, the model responded well to the seasonal changes in temperature (Figure 4, showing only monthly cohorts), with EF falling between 0.44 and 0.94, except for the Beaufort Sea. The model performed poorly in Beaufort Sea, where the growth patterns seemed more distinct between different monthly cohorts (Figure 4). The modeled lengths-at-age in Beaufort Sea were underestimated for the early hatchers with longer growth period and overestimated for the late hatchers. For the other regions, higher EF per each region in the 2005 validation simulations was correlated to different lengths-at-hatch, namely, 4 mm for Baffin Bay (EF = 0.80), 6 mm for Kitikmeot region (EF = 0.94), and 7 mm for Laptev Sea (EF = 0.76; Table S1).

Modeled growth was sensitive to the group of parameters involved in food consumption, when these were individually varied with ±20% of their nominal values (Figures 5 and S15, Text S2, Table S2). The fraction P of C_max and CD_n resulted in a proportional change in length (▵L) of ±0.08 after 30 days of growth, while CD_f induced a similar but negative response in ▵L of 0.08 and 0.05, respectively. The lower changes were induced by the temperature and CD_e with a proportional ▵L from the reference length <0.01 when the parameters values were increased (Figure 5B). When the values were lowered, however, the most sensitive parameter was temperature, resulting in negative ▵L of 0.13 of the reference length (Figure 5B).

When a global sensitivity analysis was applied to test the combined effect of the temperature and food-related parameters, the same group of parameters sensitive to the IPP induced the largest changes in model output (Figures 6A–E and S16–S18, Text S3). When both temperature and food-related parameters were varied, with values randomly drawn from a normal distribution, modeled larval length by the BEM showed a significantly negative trend across the temperature range tested (Figure 6A, Loess regression line). Among food-related parameters, P and CD_n showed a significantly positive linear effect on modeled length (Figure 6B and D). Each parameter contributed between 10% and 39% of the overall response variance (R² = 77%, metrics were normalized to sum 100%), with CD_n and P having the largest contribution (Figure 6F).

When metabolic rates of larvae were simulated over the range of temperatures found in the High Arctic and tested in experiments (up to 9°C; Koenker et al., 2018), growth rates were more sensitive to and increased sharply across the negative temperatures, followed by the gradually decreasing trend due to increase in respiration at higher temperatures (Figure 7A). At negative temperatures, growth rates were less sensitive to variations in the proportion P of the maximum daily ration (P = 0.73 and ±20% range for values of P = 0.58–0.87). Modeled growth fell within the range of laboratory estimates at higher temperatures, considering variations in the parameter P. At temperatures <3°C, our modeled growth rates were 2–3 times higher than laboratory estimates. To maintain a positive growth rate, a larva of 12-mm length requires a minimum value of P = 0.12–0.15 over the negative temperatures (Figure 7B). For the positive temperatures, the minimum value of P required for positive growth increased linearly with temperature from 0.12 at 0°C to 0.45 at 10°C.

Modeled length by the BEM for a 30-day-old larva was in the upper range of field estimates and 1.5–4.5 mm higher than the growth observed in laboratory experiments (Figure 8). When yolk reserves were present, which lasted about 6–13 days at 0°C, depending on the length-at-hatch, the model overestimated growth. This output is in comparison with previous observations of yolk-sac duration in Arctic cod larvae from hatching until complete absorption, which lasts about 20 days at 1°C (Aronovich, 1975; Graham and Hop, 1995). In the post-yolk period, variations in modeled growth and resultant length-at-age were related to both changes in food quantity (through variations of P) and quality (through variations in caloric density of prey). When a 25% lower caloric diet based on Artemia nauplii (1,295 cal g^–1; Vanhaecke et al., 1983) was simulated instead of a richer Calanus nauplii diet (1,716 cal g^–1; Thanassekos and Fortier, 2012, and references therein), which is implicitly used in the BEM, a 10% decrease in length-at-age resulted in the model output. This value remains higher than laboratory observations on larval Arctic cod growth. Reducing the daily ratio of food consumption in the model by 50% resulted in similar growth rates as laboratory estimates.

The sensitivity analysis of the model parameters indicated that those describing the feeding process are responsible for the largest changes in model predictions, a finding consistent with other bioenergetic models (Lough et al., 2005; Kristiansen et al., 2007). This finding is due, first, to order-of-magnitude higher food consumption rates compared to respiration rates and metabolic losses. Higher sensitivity in our modeled growth related primarily to variations in the proportion P of the maximum daily ration (C_max). Second, the energy of prey ingested by the larvae is not only a matter of quantity but also of food quality, reflecting the energetic conversion efficiency from prey to fish tissue (Støttrup, 2000).

In the BEM formulation, the energetic conversion efficiency was included as a ratio of the caloric content of Calanus nauplii/copepodites as prey to the caloric content of fish (Thanassekos and Fortier, 2012). At a lower energetic content of food, this ratio is reduced, hence less energy is gained. A diet of enriched Artemia nauplii, commonly used in laboratory experiments as food source for fish larvae (Graham and Hop, 1995; Koenker et al., 2018), has about 25% lower caloric content than a diet comprised of Calanus nauplii and copepodites (Vanhaecke et al., 1983; Sorgeloos et al., 2001), the dominant food source for larval Arctic cod in the High Arctic (Kohlbach et al., 2017; Bouchard and Fortier, 2020). Simulating this reduction in the caloric intake in the model resulted in a 10% decrease in the length of the early-feeding larvae at 30 days posthatching. Variation in growth related to food type has also been observed in yearling (age 1+) Arctic cod. Fish fed Calanus copepods had significantly higher growth rate compared to other food sources such as Themisto amphipods or capelin filets (Hop et al., 1997). Such diet differences may have important consequences for larval development. The consequences of poor nutrition during larval development may be obvious, for example, through lower growth and survival, deformities or malpigmentation (Shields et al., 1999; Payne and Rippingale, 2000), but in many cases, the consequences may be obscure, as in effects on temperature tolerance or limited growth during later life stages (Støttrup, 2000).

First-feeding Arctic cod larvae consume mainly copepod eggs and nauplii (80% of the stomach content) for the first days to weeks, switching gradually to a diet that includes copepodites with increased larval ability to capture larger prey size (Walkusz et al., 2011; Bouchard and Fortier, 2020). A ratio of 1/10 was suggested between larva attack speed (10 SL s^–1 estimated for Atlantic cod) and prey escape speed (100 SL s^–1) in the calculation of capture success (Fiksen and MacKenzie, 2002). This ratio simulates a prey-size preference of <5% of larval body length, in agreement with observations on Arctic cod diet by Bouchard and Fortier (2020) and on Atlantic cod diet by Munk (1997). As long as this ratio is maintained, swimming speed of both predator and prey would increase with size, resulting in little difference in the outcome of any encounter; hence, the balance between energetic spending on foraging and energetic gain on prey captured would presumably be maintained (Folkvord, 2005; Lough et al., 2005).

At high prey density, feeding success increases while the time spent foraging decreases, hence less energy spent (Puvanendran and Brown, 1999; Monk et al., 2006). We assumed an average daily feeding ration over the growth season mediated through the parameter P, which was fitted to best simulate the observed growth of wild larvae in four High Arctic regions, which represent the survivors at the end of summer. Our model output reflected well field observations, even though it relied on a fixed P value per region–year; hence, it did not account for the spatiotemporal variations in prey density. In reality, suboptimal prey concentrations during the first feeding result in slower growth, which cannot be sustained for more than 2–3 weeks, resulting in higher mortality (Monk et al., 2006; Kristiansen et al., 2007; Koenker et al., 2018). Therefore, isolated records on positive growth are not necessarily indicative of longer term consequences or survival in a patchy prey environment (Laurel et al., 2018). While it seems intuitive, larval growth is not always correlated with prey density in the natural habitat (Michaud et al., 1996). Increasing prey density, on the other hand, can lead to better survival because larvae prioritize their predator avoidance while still maintaining growth (Fiksen and Jørgensen, 2011). Therefore, predicting survival, unlike growth, needs to consider the patchiness of prey and predators in the natural environment as well as larval behavior related to predator avoidance.

Environmental temperature is the most important variable affecting metabolic rates of fish larvae (Ikeda, 2016). In theory, larval growth increases with temperature as long as prey occur at sufficient density to offset the increased metabolic costs (Clarke, 2006). In our BEM, food consumption increases asymptotically with temperature, based on estimations by Thanassekos and Fortier (2012) on the temperature effect on stomach fullness in wild fish larvae, while respiration increases linearly, reaching a maximal modeled growth at an optimum temperature of 0.7°C. Laboratory estimates, on the other hand, indicate an optimal growth temperature of 5°C (Koenker et al., 2018). The large difference in optimum temperature is not only a result of model formulation but also represents an image of data used in model parametrization. Wild larvae sampling covered only a spatiotemporal snapshot of the High Arctic, while experimental work has its own limitations, such as temperature range, duration, or acclimation potential. These aspects should be kept in mind in model evaluation and comparison. When the same temperature range will be covered by field sampling and experimental work, then different model formulations can be better tested and compared.

Modeled growth by BEM in all four High Arctic regions was sensitive to variations in temperature as small as 1°C. As a result, the model predicted different growth curves for larvae with different hatch dates. This aspect is important to consider, as the model needs to capture accurately the seasonal increase in temperature from near freezing at the beginning of the hatching period to barely above zero when simulating growth of larvae in the High Arctic. However, field observations on wild larval growth were limited to High Arctic regions, while similar observations from warmer regions remain scarce. Laboratory experiments concur with the increase growth with temperature, but their key limitation resides in the difficulty of simulating the lower thermal range observed in the natural habitat, that is, the subzero temperatures (Graham and Hop, 1995; Koenker et al., 2018). Nevertheless, laboratory experiments emphasized the need to account for both species- and stage-specific differences in the thermal response of growth models (Koenker et al., 2018). Consequently, models are useful for bridging these thermal gaps between field and laboratory observations.

While the relationship between respiration and temperature is relatively well known in fish, suggesting a positive and monotonic increase, close to exponential in shape (Clarke, 2006), estimations of food consumption with temperature remain scarce. This relationship is even more important to define, as variations in food consumption with temperature can largely dictate the optimum thermal window and growth patterns in fish larvae, as shown in our model. At higher than optimum temperatures for growth, exponential increase in respiration rates results in hyperventilation (Drost et al., 2014, 2016), while food uptake is increasingly used to cover the exponentially higher metabolic demands, finally resulting in declining growth rates (Brett, 1979).

Surprisingly, observations on immature adult Arctic cod (age 1+ and 2+) indicated that growth rate did not vary significantly over the temperature range 0°C–8°C, while respiration was strongly affected by temperature (Kunz et al., 2016). This finding implies that fish could have indeed compensated for increased respiration with higher food uptake. This assumption fits our results on the minimum proportion P of the maximum daily ration required to maintain positive growth. Nevertheless, this strategy for growth would then suggest that at least immature adult Arctic cod have the physiological potential to expand their distribution southward into warmer areas, given that enough food is available in the environment to compensate for increasing metabolism related to higher temperature. Although when the narrow range of thermal tolerance of Arctic cod juveniles overlapped with that of subpolar fish (e.g., saffron cod Eleginus gracilis), growth and lipid accumulation exceeded that of their southern congeners only at water temperatures <2°C (Laurel et al., 2016; Copeman et al., 2020), indicative of limiting factors of Arctic cod distribution other than temperature and food availability (probably food type and quality). This compensatory strategy for growth would further imply that at low temperature, Arctic cod require less food than at higher temperature or than southern competitor species to maintain similar growth rates, thus presenting an advantage in the High Arctic (Laurel et al., 2018).

Generally, fish larvae from cold ecosystems are believed to require less food than larvae from warm ecosystems (Pörtner et al., 2005). This strategy does not concur with our modeled growth in four Arctic regions, which showed that larger daily ration (mediated through parameter P) was needed in regions with lower temperature to reproduce the field-observed growth. A limitation in our models arises from the data used for parameter fitting and validation, which was restricted to the High Arctic region, north of the Arctic Circle, with a surface temperature range from –1.8°C to 5°C at best (Figures 1 and S1–S4). Furthermore, larvae are more thermally sensitive than juvenile and adult Arctic cod (Drost et al., 2016; Laurel et al., 2017) and could have a different strategy for optimal growth or for survival, upon which we can only speculate. Clear, however, is that these larvae have an adaptation potential that makes them thrive at near-freezing temperature.

Ultimately, the combined effects of increasing metabolic demands with temperature and compensation with food consumption are not only a matter of food quantity but of energetic gain. At low temperature, fish have lower metabolic demands, spend less energy foraging but have high energetic gain when a diet of lipid-rich Arctic prey is consumed (Kohlbach et al., 2016, 2017). Therefore, even when PA is limited, low temperatures could still yield the highest growth for Arctic cod in the High Arctic. Our model indicated that growth was less sensitive to variations in P at negative temperatures, implying that a relatively constant growth could be maintained despite spatiotemporal variations in PA, likely because at low temperatures, slow digestion limits food intake (Hop and Tonn, 1998). At high levels of available prey, maximum growth is normally expected to occur at higher temperatures (Lough et al., 2005). This dynamic response of metabolic processes to temperature agrees with laboratory and field observations on the growth of other fish species (Brett and Groves, 1979, Otterlei et al., 1999). For Arctic cod larvae, however, such growth dynamics remain speculative because of the lack of direct observations.

Modeled growth rates in larval Arctic cod resulted in similar rates to those observed in the wild and estimated from otolith microstructure (Bouchard and Fortier, 2011; Bouchard et al., 2017), except for the yolk-sac period when growth was overestimated by the model, predicting a shorter duration (Graham and Hop, 1995; Kent et al., 2016). Of interest is that the experimental work published to date on the growth of larval Arctic cod reported lower growth rates than those estimated by our model or observed in the wild at similar temperatures (Graham and Hop, 1995; Koenker et al., 2018; Bender et al., 2021). Our BEM might overestimate larval growth compared to laboratory settings for two reasons. First, food quality is better in nature (by 25% higher caloric content), with a variety of lipid-rich prey being available in the High Arctic (Kohlbach et al., 2016). Second, the model was adjusted using data from wild age-0 fish collected between August and October; that is, the survivors of the larval cohorts remaining after selective mortality likely removed the smaller-at-hatch and slower growing individuals from the population through predation and starvation (Peck and Hufnagl, 2012; Jørgensen et al., 2016). Selection for large size-at-hatch has been documented in larval Arctic cod (Pepin et al., 2014), while selection for fast growth is likely to occur among larvae hatching at the same time. Faster growth in wild larvae compared to reared larvae has been observed in other species, notably Atlantic cod (Finn et al., 2002; Folkvord, 2005).

To date, we are unaware of bioenergetic models of larval Arctic cod other than the one presented by Thanassekos and Fortier (2012), which derived daily ration estimates from immature adults in the absence of data on larvae. In the absence of accurate field prey data, model approximations can be informative despite their limitations, in the same way as experimental work that uses two fixed settings of prey density (low vs. high scenarios) over the growth period (Koenker et al., 2018). Given the variability in model validation studies and the uncertainty in model inputs (such as oxycalorific equivalents and the energy content of predators and their prey, among others), testing models at all empirically defensible thermal, size, and feeding levels, even in the laboratory, may not be possible. Our models, based on average daily values of P and regional temperatures, might have masked the potential for patchiness in prey concentration and temperature that could potentially boost Arctic cod consumption and growth or, on the contrary, lead to starvation.

Feeding behavior continues to be a difficult process to model and questions related to prey preference and foraging are particularly challenging (Fiksen and MacKenzie, 2002). Equations including individual foraging behavior (Kristiansen et al., 2007; Peck and Daewel, 2007) and depth-varying conditions (Fiksen and MacKenzie, 2002) or coupled biophysical simulations on match–mismatch feeding dynamics (Kristiansen et al., 2009) would perhaps account for small-scale variability in prey abundance and temperature. The increased accuracy in predicting growth provided by more sophisticated schemes, however, would come at higher computational cost. Moreover, field observations of prey data are required for such simulations, which are largely unavailable in the High Arctic, especially during the ice-covered period, therefore limiting the capacity of such a model to resolve, for example, the impact of hatch dates on the growth trajectory of different cohorts at a circum-Arctic scale. Alternatively, prey data could indeed be estimated from biochemical NPZD (nutrient, phytoplankton, zooplankton, and detritus) or ecosystem models, which come, however, with additional prediction errors that need to be considered. Thus, the key advantage in simplifying the feeding processes in a bioenergetic model, such as the BEM presented here, relies on the possibility to be further applied at larger spatial scales to test larval growth dynamics with limited input data.

The bioenergetic model presented here provides a robust description of the general pattern in larval Arctic cod growth, synthesizing metabolic processes that are of both theoretical and practical importance to modelers and physiologists. With limited information on feeding behavior and underlying relationships between metabolic processes and temperature, BEM proved useful in estimating growth of larval Arctic cod, well within the range of field observations.

Modeling the variability of individual growth and its determining factors can considerably increase our understanding of Arctic cod recruitment, population dynamics, and response to environmental changes. Our model allowed us to test hypotheses with only limited information at hand on prey availability or feeding dynamics and to identify crucial knowledge gaps, such as the need for estimating daily ration with larval weight and temperature and yolk-consumption dynamics. Further consideration should also be given to different types of food used in experimental work. Aquaculture feed typically is lower in quality than Arctic copepods and may lead to suboptimal growth even with high availability. Ultimately, using modeling tools to bridge gaps between laboratory and field studies can prove crucial in expanding our understanding of important species in Arctic seas.

The R scripts generated in this study for model formulation, data analysis, and figure plotting are publicly available on Github and can be accessed at https://github.com/carmenldavid/ArcticCod_BioenergeticModel.

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

Figure S1. Sea surface temperature (SST) in Baffin Bay.

Figure S2. Sea surface temperature (SST) in Beaufort Sea.

Figure S3. Sea surface temperature (SST) in Kitikmeot region.

Figure S4. Sea surface temperature (SST) in Laptev Sea.

Text S1. Fitting parameter P and model validation.

Figure S5. Fitting parameter P and model validation: Baffin Bay 2006.

Figure S6. Fitting parameter P and model validation: Baffin Bay 2008.

Figure S7. Fitting parameter P and model validation: Beaufort Sea 2010.

Figure S8. Fitting parameter P and model validation: Beaufort Sea 2011.

Figure S9. Fitting parameter P and model validation: Beaufort Sea 2014.

Figure S10. Fitting parameter P and model validation: Kitikmeot 2006.

Figure S11. Fitting parameter P and model validation: Kitikmeot 2011.

Figure S12. Fitting parameter P and model validation: Kitikmeot 2015.

Figure S13. Fitting parameter P and model validation: Laptev Sea 2003.

Figure S14. Fitting parameter P and model validation: Laptev Sea 2007.

Table S1. Model efficiency in predicting length-at-age in 2005 for different larval sizes at hatch.

Text S2. Sensitivity of growth to individual model parameters.

Table S2. Variation in model parameters and contribution to prediction error (R²) shown by two larval sizes.

Figure S15. Sensitivity analysis by individual parameter perturbation (IPP).

Text S3. Sensitivity of growth to the combined effect of model parameters.

Figure S16. Pearson correlation between modeled larval length (small larvae of 6 mm) and model parameters.

Figure S17. Pearson correlation between modeled larval length (large larvae of 12 mm) and model parameters.

Figure S18. Parameters ranking from a global sensitivity analysis.

This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund, through the Ocean Frontier Institute.

Authors declare no competing interests.

• Contributed to conception and design: CLD, RJ, JH.

• Contributed to acquisition of data: CLD, CB.

• Contributed to analysis and interpretation of data: CLD, CB.

• Drafted and/or revised this article: CLD, RJ, CB, HH, JH.

• Approved the submitted version for publication: CLD, RJ, CB, HH, JH.

How to cite this article: David, CL, Ji, R, Bouchard, C, Hop, H, Hutchings, JA. 2022. The interactive effects of temperature and food consumption on growth of larval Arctic cod (Boreogadus saida): A bioenergetic model. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2021.00045

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

Guest Editor: Dominique Robert, Institut des Sciences de la Mer, Université du Québec à Rimouski, Canada

Knowledge Domain: Ocean Science

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