This study aims to understand how ocean warming caused by climate change could alter egg survival and limit growth potential in cod species (Polar cod, Atlantic cod, and Greenland cod) found off the east coast of Canada. The analysis focused on the projected changes in temperature in ocean shelf areas where ocean depth is ≤400 m from 40°N to 70°N and 20°W to 70°W. Grid cells off the coast of Greenland were excluded. Climate influence on water temperature was projected for near-surface (∼5 m) and bottom shelf water (≤400 m) to account for gadids with demersal eggs (Greenland cod) and those with pelagic eggs (Atlantic and Polar cod).

To investigate the projected impacts of climate change on water temperature, historical baseline temperatures were calculated using the European Center for Medium-Range Weather Forecasting (ECMWF) Ocean Reanalysis Pilot 5 (ORAP5; Zuo et al., 2015). Monthly averages of bottom and surface seawater temperatures were calculated for a 25 years historical baseline period (1981–2005). Climate hindcasts and projections were then derived from simulations of CMIP5 global climate models (GCMs) with a single ensemble run (r1i1p1) selected for each model (see Supplementary Table 1 for model information). This analysis focuses on the RCP 8.5 greenhouse gas emissions scenario, which is the highest of the RCP scenarios and can be interpreted as “business as usual”. Any projected changes on our results from a more moderate emissions scenario (e.g., RCP 4.5) are typically found to be similar in trend but smaller in magnitude than RCP 8.5. Before the analysis, each GCM was regridded using nearest neighbor interpolation to match the ORAP5 grid resolution (1° × 1°). GCMs contain systematic errors (biases) in their output due to coarse resolution, simplified physics or thermodynamics, or incomplete understanding of the climate system. To overcome this problem, climate models can be bias corrected using various approaches. Here, we apply the delta-change approach (Hawkins et al., 2013) for bias correcting the GCMs (see Supplementary Figures 1, 2 for GCM ensemble temperature predictions and debiased GCM ensemble predictions over the study period). For each GCM, monthly averages for the historical baseline period (1981–2005) were calculated and subtracted from the monthly temperatures for each year (1981–2100) to create deltas. These differences, or deltas, were then added on to the observed baseline for each month (ORAP5; monthly averages from 1981 to 2005). Due to the non-linear response of growth and survival across temperature, we applied the thermal developmental models (i.e., temperature vs. egg survival and growth potential; Table 1; Supplementary Figure 3) to relevant monthly water temperature projections before aggregating values. Near-surface temperatures were applied to all egg survival models, with the exception of Greenland cod, for which bottom temperatures were used for its demersal eggs. Since shallow habitats are important rearing areas for all three species, near-surface temperatures were used exclusively for growth models. Each development metric was then averaged over baseline (1981–2005), near future (2026–2050), mid-century (2051–2075), and end-of-century (2076–2100) climatology periods. The CMIP5 ensemble median was taken across all models to generate growth potential and egg survival surfaces for the study region. To assess uncertainty of the projections, the temporal change between the baseline period and each simulated future period (Δ egg survival/growth potential) was compared to the ensemble spread within a simulated time period (as in Dahlke et al., 2018). The ensemble spread was calculated as the median absolute deviation (MAD) across ensemble members for each time period. Model results were not considered robust where the temporal change in a metric was less than the spread of that metric within a time period (i.e., Δ egg survival/growth potential/MAD <1 as in Dahlke et al., 2018). Egg survival was normalized (from 0 to 1) to compare the relative changes in survival across species with changes in temperature.

Self-sustaining populations require seasonal environmental conditions that support multiple developmental life stages. Toward a more representative index of regional population health, we combined the projected surfaces of both egg survival and juvenile growth potential into a single habitat suitability index (Mazzotti et al., 2007). Specifically, the projected surfaces for each development metric and species were normalized (from 0 to 1) by dividing each grid cell by the maximum potential value of the development model. The geometric mean of egg survival and growth potential were calculated for each cell across time periods to get an overall habitat suitability value from 0 to 1 with values of 0 being not suitable for egg survival or growth and values of 1 being highly suitable for both metrics.

To compare the habitat suitability surface to the known distributions of our study species, distribution polygons were created by aggregating point occurrence data for Polar cod (GBIF.org, 2021a), Atlantic cod (GBIF.org, 2021b), and Greenland cod (GBIF.org, 2021c) from Global Diversity Information Facility (GBIF). These aggregations were performed using the Aggregate Points tool in ArcGIS Pro (V2.5.0) and the aggregation tolerance parameter was set to 1 decimal degree. These aggregations were then smoothed using the Polynomial Approximation with Exponential Kernel (PAEK) smoothing algorithm and the tolerance set to 2 decimal degrees. Areas overlapping land were clipped from the smoothed polygons.

We summarized climate-related shifts in latitude using distribution centroids. The distribution centroid gives the latitude at which the sum of grid-cells below and above that point are equal and as a result can be used to detect spatial shifts in the distribution of values (see e.g., Cheung et al., 2009). The distribution centroid (C) calculated as:

where w_i and Lat_i are the value and the mean latitude at each cell i and N is the total number of cells. Values were determined for the egg survival, growth potential, and the combined thermal habitat suitability surfaces across all time periods and species. For all metrics and species, the distribution centroid for each projected time period was subtracted from the historic baseline period to quantify potential distributional shifts with warming.

To quantify changes in areas of high thermal suitability (HTS), the area of grid cells above the 50^th percentile of baseline values was calculated for each metric over time. The percentile threshold was arbitrary but provided a useful benchmark with which to compare change. Areas for each grid cell were approximated using the area function in the raster package (Hijmans, 2020) in R. This function approximates the area by multiplying the height of the cell by the width of the cell at the latitudinal middle of the cell. The data was then arranged in order of increasing cell values for each species and metric. The 50^th percentile was taken as the cell value at which the cumulative proportion of the area was 0.50. The sum of the area above this threshold was calculated across all climatological periods. Analyses were carried out using Python (V.3.6.9), R (V.4.0.2), and ArcGIS Pro (V.2.5.0).

Temperature has been widely used to predict species distribution changes associated with climate change (e.g., Morley et al., 2018; McHenry et al., 2019). However, mechanism-based risk assessments targeting vulnerable life stages in a climate-scenario context are rare, especially for marine species inhabiting Arctic regions. Dahlke et al. (2018) forecasted changes in Atlantic cod and Polar cod distributions in the northeast Atlantic based on shifting embryonic thermal tolerance to ocean acidification. This study builds upon Dahlke et al. (2018) by expanding consideration of seasonal conditions to support two critical early life history stages of gadid species in the northwest Atlantic. This study also follows on climate vulnerability assessments that emphasize the need to consider multiple life-history constraints in forecasting fish distributions under changing thermal environments (e.g., Hare et al., 2016).

Gadids are behaviorally responsive to temperature and have historically altered distributions according to climate conditions. Atlantic cod, the most well studied of our focal species, has shifted distributions poleward with warming in both the northeast Atlantic (Drinkwater, 2006; Engelhard et al., 2014) and northwest Atlantic (Rose et al., 1994). Southern shifts have also been documented during cold water anomalies in the late 1980s and early 1990s in the northwest Atlantic (Rose et al., 2000). Historical spatial shifts were extensive and rapid in some cases, as exhibited by the northern incursions of commercially viable populations off west Greenland during the warm period of the 1920s and 1930s (Drinkwater, 2006) and the collapse of the Labrador fishery as the northern (Atlantic) cod stock moved to the south (Rose et al., 2000). Distributional shifts of Greenland cod in the northwest Atlantic have not been documented, but populations in the Pacific have made extensive shifts poleward into the North Bering Sea following recent warming and loss of sea ice (Stevenson and Lauth, 2019). Further south in the Pacific, a series of marine heatwaves beginning in 2014 led to negative impacts on G. macrocephalus populations that resulted in a complete closure of the fishery in 2020 (Barbeaux et al., 2020; Laurel and Rogers, 2020). Populations at the historical southern range of this species in the northeast Pacific may be extirpated (Palsson, 1990). Less is known about environmental-driven distribution shifts of Polar cod, but they were also observed to shift north in response to the warmer waters of the 1920s and 1930s (Drinkwater, 2006) and have also shifted distributions to southern areas during cold water conditions in the northwest Atlantic (Drinkwater and Mountain, 1997; Rose et al., 2000). Seasonal distributions of Polar cod also change with the presence of cold water (Mueter et al., 2016; Baker, 2021), and models in the Barents Sea suggest a northeast spawning assemblage shift with warming (Huserbråten et al., 2019).

The effects of climate change are manifested in complex ways across species and life stages, based on life history, local geography and oceanography, and interspecific ecological interactions. Not surprisingly, we project suitable habitat for all considered life stages and species to shift poleward by 2100 with non-uniform distribution shifts across species (Fulton, 2011; Hazen et al., 2012; McHenry et al., 2019). However, the changes to available habitat area within these ranges varies broadly among life stages. For example, the area of habitats most conducive to egg survival of Polar cod are strongly reduced over time (55% of baseline area), whereas they are expected to expand for Atlantic cod and Greenland cod (153 and 127%, respectively). In contrast, the area of habitats with strong growth potential will decrease slightly for Polar cod (72% of baseline area) but strongly decrease for the other two species (Atlantic cod: 26% of baseline area; Greenland cod: 37% of baseline area). The relatively less dire growth predictions for Polar cod are supported by Bouchard et al. (2017) who predict growth conditions for juveniles to improve initially as climates warm through 2100. The negative projections for egg survival of Polar cod are also consistent with models applied to the northeast Atlantic (Dahlke et al., 2018), in that the northward shift of the southern range limit means Polar cod will be progressively squeezed within a smaller area of suitable habitat. In contrast, boreal species like Atlantic cod and Greenland cod can shift distribution limits northward, encroaching into areas previously occupied by Polar cod (Arctic Council, 2010; Stortini et al., 2016; Huserbråten et al., 2019). However, unlike the northeast Atlantic, the areal extent of growth habitats for the boreal/subarctic species examined in this study will be reduced as fish move beyond latitudes with expansive shelf areas (e.g., Grand Banks, Gulf of St. Lawrence; Morley et al., 2018).

Beyond geography, our study species vary in their vulnerability to warming. Polar cod, for example, have relatively narrow thermal tolerances (Supplementary Figure 3), extended pelagic larval phases (Walkusz et al., 2011), and strong associations with sea ice (Bouchard and Fortier, 2011; Huserbråten et al., 2019). These characteristics increase the likelihood of encountering unfavorable environmental conditions, particularly given that sea ice extent is declining in the Arctic (Stroeve et al., 2012) and is projected to be nearly ice free in the summer by mid-century (Wang and Overland, 2012; Notz and SIMIP Community, 2020). However, forecasted climate velocities at subsurface layers are higher than surface layers (Brito-Morales et al., 2020), which will be especially challenging for life stages of species with narrow thermal tolerance like the eggs of Greenland cod. The broader thermal tolerance and extended spawning period of Atlantic cod (Table 1; Ings et al., 2008) also provide some safeguards against environmental anomalies compared with other gadid species.

The poleward shifts in spawning habitat for gadids in the western Atlantic (1.0°–3.2° latitude) in this study are not as dramatic as those predicted for Atlantic cod and Polar cod in the eastern Atlantic (Dahlke et al., 2018). For Atlantic cod, our models predict relatively high spawning suitability in their current spawning areas off Newfoundland (50° latitude) whereas models in the northeast Atlantic predict a near absence of spawning for Atlantic cod south of the Arctic Circle by 2100 under RCP 8.5 scenarios (Dahlke et al., 2018). The different projections in the northwest Atlantic are likely driven by the region’s distinct oceanography that includes transport of much colder Arctic water southward to Newfoundland. Smaller distribution shifts were also predicted for fish species of the Canadian Atlantic coast compared with other parts of North America (Morley et al., 2018). Warming oceans, associated with multi-decadal climate cycles, have previously caused large changes in marine ecosystems in the north Atlantic and such incidents can be instructive to projecting effects of future climate change. For example, an ocean warming phenomena in the 1920s–1930s in the North Atlantic corresponded with range shifts, new spawning areas, and large increases in abundance of Atlantic cod in west Greenland in addition to a host of other ecosystem changes (Drinkwater, 2006). The origin of the invading Atlantic cod was thought to be from the east (Iceland) instead of the south (Labrador) where temperature is moderated by the strength of the cold Labrador current and evidence of the broad ecosystem change was weaker compared with the northeast Atlantic (Drinkwater, 2006). Given the prevailing currents that can transport fish and larvae from southern Greenland into Baffin Bay of the Canadian Arctic (Coté et al., 2019), it is possible that Atlantic cod will emerge in the Canadian Arctic from the northeast rather than the northwest Atlantic under climate warming scenarios.

Past warming events also show that many ecosystem changes accompany warming water, including changes to primary and secondary productivity, altered species assemblages, shifts in distribution (e.g., Carscadden and Nakashima, 1997) and migration timing, and concomitant effects on fisheries and economies (Drinkwater, 2006). Our predictions, while focused on early life stages that are likely most vulnerable to direct temperature effects, do not factor additional stressors, e.g., changing salinity (Spencer et al., 2020), ocean acidification, altered harvesting pressure (Opdal, 2010; Engelhard et al., 2014), changing prey availability (Rose et al., 2000), competitive interactions (Linehan et al., 2001), or low dissolved oxygen (Pörtner and Farrell, 2008). For example, sea ice decline associated with climate change will modify ocean salinity (Li and Fedorov, 2021) and could impact egg and larvae development of gadids (Nissling and Westin, 1997; Sakurai et al., 1998; Spencer et al., 2020). Other ecosystem components or processes may contribute directly or indirectly to successful spawning and reproductive output (Bouchard et al., 2017; McHenry et al., 2019) such as the seasonal light cues required for synchronizing spawning (Van Der Meeren and Ivannikov, 2008), sea ice for egg retention and ice-algae production (Bouchard and Fortier, 2008), or other environmental covariates required to match primary and secondary prey production with first-feeding success (Laurel et al., 2021).

Climate-induced changes to prey and predator fields could also have important impacts on young life stages (Donovan et al., 2017). Within our study species, competition would most likely occur between Greenland cod and Atlantic cod that currently occupy similar habitats as juveniles (e.g., Coté et al., 2013) and prey on one another where their ranges overlap (Linehan et al., 2001). Although Atlantic cod were noted to have improved growth efficiency over Polar cod under warming scenarios (Kunz et al., 2016), competition for food between Polar cod and Atlantic cod was considered to be minimal in Svalbard where both species are found (Renaud et al., 2012). Predation may be more likely to inhibit Polar cod where they co-occur with our other two species (Renaud et al., 2012; Bouchard et al., 2017) as both Atlantic cod and Greenland cod will predate Polar cod (Scott and Scott, 1988). While spatial and temporal shifts in spawning and rearing habitat are becoming increasingly documented in cod populations (e.g., Rose et al., 2000; Drinkwater, 2006; Rogers et al., 2019), the consequences of these changes on the population are not fully understood. Our supplementary analyses (Supplementary Figure 4) indicate that the average optimum month for spawning temperature, within the range of months used in our modeling, will shift from winter to early spring over our study period. In contrast, optimum growth periods, that currently occur in late summer, remain relatively static over the study period.

Distributional changes in the gadid assemblage also introduce a number of potential impacts on other components of the ecosystem. At the juvenile stage, gadids are important forage species for fish, invertebrates, seabirds, and marine mammals (Scott and Scott, 1988; Welch et al., 1993; Linehan et al., 2001; Benoit et al., 2010). Since nutritional value (e.g., lipid energy) varies significantly among species, ontogenetic stage (Copeman et al., 2020), and thermal habitat (Copeman et al., 2017; von Biela et al., 2019), substitution of an Arctic species with a boreal species may not preclude ecosystem impacts. For example, Polar cod are the most lipid rich among these gadids at the juvenile stage (Copeman et al., 2020) and serve an energy-efficient link between secondary production and higher trophic levels. In the Beaufort Sea, models suggest energetic contributions from Polar cod would have to be replaced by other energy rich sub-Arctic species such as capelin (Hoover et al., 2013; Walkusz et al., 2013).

Establishment of gadid populations in new areas will require a consistent supply of juveniles and sufficient food, in addition to suitable environmental conditions that can vary by life stage (Rindorf and Lewy, 2006). The spatio-temporal overlap of suitable spawning and juvenile habitat is an overlooked but important component of climate resilience. In particular, new habitats must sustain embryonic, larval, and juvenile stages or at least allow for dispersal toward adequate nursery areas (deYoung and Rose, 1993; Bradbury et al., 2008), as survival during these early lifestages likely determines year class strength and recruitment to the adult population (Houde, 2008). It is reasonable to assume that population productivity will be negatively impacted if favorable habitats for different life stages are further mismatched in time and space by climate forcing.

Our habitat projections rely on several key assumptions. For example, an exogenous feeding juvenile must have sufficient food to achieve its growth potential for a given temperature. In cases where food is limited, optimal growth temperatures are typically cooler (Jobling, 1995). Indeed, where temperature-dependent growth has been investigated across species, average growth along temperature clines tends to be weaker than predicted by metabolic theory, largely due to local ecosystem dynamics (van Denderen et al., 2020). Thus it is likely that our projections will overestimate juvenile habitat suitability in some areas. In contrast, eggs do not rely on external food sources and newly hatched larvae have yolk reserves to buffer times of low productivity. We therefore have higher confidence in habitat predictions based on egg survival than juvenile growth. However, neither model accounts for existing or future species adaptations that could further influence projected habitat suitability under varying climate scenarios (Peck et al., 2018). For example, the single species-specific models we apply may not capture existing local adaptations to environmental conditions (Bradbury et al., 2010; Berg et al., 2015; Wilson et al., 2020; Spies et al., 2021) or the responses of these populations to new climate driven selection pressures. Moreover, behavioral adaptations, such as altered depth use (Kleisner et al., 2016; Poloczanska et al., 2016) or spawning chronology could also mitigate the need for latitudinal distribution shifts in response to pervasive warming. Continued advances in species distribution modeling techniques such as hybrid models (Rodríguez et al., 2019) and adaptive models (Bush et al., 2016; Chen et al., 2020) will allow for more in-depth and holistic assessments of species vulnerabilities to climate change (Peck et al., 2018).

The climate projections we incorporated to support the physiological models were effective at capturing surface temperature when compared with observations (Supplementary Figure 1), although it was essential to bias-correct the projections to obtain realistic historical ranges. GCMs have evolved considerably over the last few decades and are now able to reproduce, with high confidence, seasonal cycles and mean temperature variance on interannual and centennial cycles (Flato et al., 2013). Projected sea surface temperatures for the northwest Atlantic are closely linked with the ability of climate models to adequately simulate the extent and evolution of sea-ice and variability can have strong impacts on temperature distributions. Model evaluations suggest that climate models are able to reproduce the Arctic sea-ice cycle and sea-ice extent with an error margin of less than 10% across models (Flato et al., 2013); however, there is a wide spread in the mean state and sensitivity of warming to CO₂ emissions that introduces significant uncertainty across models (Notz and SIMIP Community, 2020). This suggests there is considerable room for improvements to the climate projections within our sea-ice dominated study region and such improvements would have an impact on the habitat suitability estimates presented here. Still, the GCM and Earth System Models (ESMs) simulations over the last century are consistent with the historical observations of global average surface temperature, ocean heat budget, and carbon cycle responses, suggesting that these kinds of models provide valuable information on the long-term large-scale climate responses (Flato et al., 2013).

Since GCMs apply simplified physics and thermodynamics at relatively coarse scales, small-scale processes such as eddies and turbulence are parameterized, whereas tides are not included and bathymetry and freshwater are represented at a coarse scale. This can lead to systematic biases in these models if not corrected with methods such as the delta-change approach used in this study. This adjustment corrects some of the bias at the regional and local level and is a critical step in using climate model output, particularly for threshold studies such as our habitat analysis. Still, the relatively coarse scale CMIP5 models do not resolve inner-coastal shelf processes such as tidal mixing, freshwater input, and bathymetry very well, resulting in higher uncertainty of the projections closer to land. Capturing local-scale dynamics would require either a dynamical or statistical downscaling and was beyond the scope of this study. Despite these uncertainties, our models largely correspond to existing species distributions and also illustrate the complex responses of species assemblages and life stages within a species to climate change.