Forecasted Shifts in Thermal Habitat for Cod Species in the Northwest Atlantic and Eastern Canadian Arctic

Cote, David; Konecny, Cassandra A.; Seiden, Jennica; Hauser, Tristan; Kristiansen, Trond; Laurel, Ben J.

doi:10.3389/fmars.2021.764072

ORIGINAL RESEARCH article

Front. Mar. Sci., 29 November 2021

Sec. Marine Ecosystem Ecology

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.764072

This article is part of the Research Topic Toward an Integrated Understanding of Ecological Responses to Variable Climate Changes and Marine Environments View all 9 articles

Forecasted Shifts in Thermal Habitat for Cod Species in the Northwest Atlantic and Eastern Canadian Arctic

$\r\nDavid Cote*$ David Cote^1*

Cassandra A. Konecny¹

Jennica Seiden¹

Tristan Hauser²

Trond Kristiansen^3,4

Ben J. Laurel⁵

¹Fisheries and Oceans Canada, Northwest Atlantic Fisheries Centre, St. John’s, NL, Canada
²Wood E&IS GmbH, Frankfurt am Main, Germany
³Farallon Institute for Advanced Ecosystem Research, Petaluma, CA, United States
⁴Norwegian Institute for Water Research, Oslo, Norway
⁵Fisheries Behavioral Ecology Program, Alaska Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Newport, OR, United States

Climate change will alter ecosystems and impose hardships on marine resource users as fish assemblages redistribute to habitats that meet their physiological requirements. Marine gadids represent some of the most ecologically and socio-economically important species in the North Atlantic, but face an uncertain future in the wake of rising ocean temperatures. We applied CMIP5 ocean temperature projections to egg survival and juvenile growth models of three northwest Atlantic coastal species of gadids (Atlantic cod, Polar cod, and Greenland cod), each with different thermal affinities and life histories. We illustrate how physiologically based species distribution models (SDMs) can be used to predict habitat distribution shifts and compare vulnerabilities of species and life stages with changing ocean conditions. We also derived an integrated habitat suitability index from the combined surfaces of each metric to predict areas and periods where thermal conditions were suitable for both life stages. Suitable thermal habitat shifted poleward for the juvenile life stages of all three species, but the area remaining differed across species and life stages through time. Arctic specialists like Polar cod are predicted to experience reductions in suitable juvenile habitat based on metrics of egg survival and growth potential. In contrast, habitat loss in boreal and subarctic species like Atlantic cod and Greenland cod may be dampened due to increases in suitable egg survival habitats as suitable juvenile growth potential habitats decrease. These results emphasize the need for mechanistic SDMs that can account for the combined effects of changing seasonal thermal requirements under varying climate change scenarios.

Introduction

Climate change is having widespread impacts on marine species (Pörtner et al., 2014; Lemmen et al., 2016) through altered ocean temperatures, sea-ice loss, decreased surface salinity, ocean acidification, and reduced oxygen availability (Stortini et al., 2016). These environmental changes affect organisms through physiological and ecological pathways (Bryndum-Buchholz et al., 2020a,b) and can result in genetic, demographic, and/or behavioral responses. Often the combined effect is observed through expansion or contraction of geographic range or through modified depth use, altered abundance, and/or changes in community structure and functioning (Bryndum-Buchholz et al., 2020a,b). These changes in turn impose hardships on commercial and cultural fisheries established to target particular fish assemblages that are no longer available to harvest (Ojea et al., 2020).

Predicting how climate change may alter the abundance and distribution of marine species is an important application of species distribution models (SDMs). The approach is rapidly expanding in the field of marine ecology (Freer et al., 2017; Robinson et al., 2017; Melo-Merino et al., 2020) and is increasingly being used to guide adaptive fisheries management, marine conservation plans, impact assessments, and policy decisions (Freer et al., 2017; Tittensor et al., 2019; Bryndum-Buchholz et al., 2020a). Species distribution forecast models combine known biological-environmental relationships with downscaled oceanographic models to project species range shifts under one or more Representative Concentration Pathways (RCPs) of greenhouse gas concentration trajectories. Broadly speaking, SDMs can be categorized as: (i) correlative models that use occurrence or abundance data and a set of environmental variables at known locations; or (ii) mechanistic models that apply an understanding of fishes’ physiological limits and biophysics. Correlative models are the most widespread in the literature and are often based on observed environmental relationships of adults from independent fish surveys conducted in a particular management area. Correlative approaches rely on multiple years of observed environmental variability, but have the disadvantage that they are based on a relatively narrow range of observed environmental conditions within the physiological optima of the species (Kearney et al., 2010). Correlative SDMs may therefore overestimate the suitability of future warming scenarios and miss critical thresholds, especially for species that have not currently experienced thermal stress in their region (Anderson, 2013). Moreover, correlative models may perform poorly when the predictive variables do not represent the mechanism limiting the species distribution (Peck et al., 2018) or are obscured by survey biases (e.g., Araujo and Guisan, 2006). In contrast, mechanistic SDMs provide an attractive alternative as they better capture a broader range of responses of species to future environments (e.g., non-linear, threshold, changes in environmental variability, etc.) and allow users to examine constraints of other lifestages and habitats that are unavailable in typical time-series survey data. Such approaches have gained traction because they can be relatively simplistic once the biological rate (e.g., growth or survival potential) is parameterized across a broad range of environmental conditions. While biological information can be limiting in rare species, such approaches have been used to forecast habitat conditions for well-studied species like gadids and other commercially and ecologically important fish (e.g., Dahlke et al., 2018).

Mechanistic SDMs generally focus on temperature given its pronounced effect on a fish’s life, particularly for early life stages that are less able to physiologically and behaviorally mitigate unfavorable conditions (Pankhurst and Munday, 2011; Dahlke et al., 2020). This is particularly true for high latitude species, which are characterized by narrower thermal tolerance ranges (Sunday et al., 2019). Thermal responses of fish are generally dome-shaped, but are unique to each species in terms of its width (thermal breadth), peak (thermal optimum), and placement along a thermal range. Moreover, there can be different thermal optima for maximum and efficient development/growth rate, which further varies by life stage. Two potential physiological mechanisms driving climate-related change in fish communities are egg sensitivity to temperature and the effects of temperature on juvenile growth. Egg stages generally have narrow temperature tolerance compared to older life stages (Dahlke et al., 2020) and deviations from thermal optima can negatively impact fertilization (Klein et al., 2017), egg size (Klein et al., 2017), hatching success (Laurel and Rogers, 2020), and likelihood of complete development of cardiovascular and other homeostatic structures (Pörtner et al., 2017). Moreover, temperature changes can go on to affect a suite of other developmental and growth processes linked to timing and magnitude of food availability (Kristiansen et al., 2011; Durant et al., 2019; Laurel et al., 2021). Juvenile fish must also budget temperature-dependent energetic demands for storage (fat) and somatic growth (Sogard and Olla, 2000). In higher latitudes, juvenile growth and size-at-age are key metrics regulating population dynamics, largely by way of early life survival (i.e., first year), maturation onset, fecundity, and offspring quality (Campana, 1996; Sogard, 1997; Swain et al., 2007). When environmental conditions are unfavorable, fish cohorts have poor survival, often succumbing to starvation and/or high levels of predation, and ultimately affecting future recruitment to adult populations (Houde, 2008; Laurel et al., 2017). The make-up of fish communities in a particular area is thus a reflection of the suitability of that area’s environment for various fish species.

The northwest Atlantic Ocean comprises vast productive fishing grounds but is experiencing regional climate-induced changes to its oceanographic conditions with some regions forecasted to rise almost three times faster than the global average (Pershing et al., 2015; Saba et al., 2016). Within inshore zones, gadids such as Polar cod (Boreogadus saida), Atlantic cod (Gadus morhua), and Greenland cod (Gadus macrocephalus) are prominent species in terms of their ecological and/or socio-economic value. Atlantic cod is the most southerly distributed of the three species and was/is the dominant fish predator in the northwest Atlantic, and currently supports fisheries worth $20–30 million (historically over 100 million) (Fisheries and Land Resources, 2019). Polar cod are distributed furthest north and are an important food source for a variety of Arctic fish, marine mammals, and seabirds (Welch et al., 1993; Harter et al., 2013). Intermediate in latitudinal distribution to the aforementioned species is Greenland cod, which although not commercially fished in Canadian waters, is an important subsistence food source for Indigenous communities on the coast of Labrador (Dombrowski et al., 2013). Recent research suggests that warming ocean temperatures may greatly reduce the available habitat of Polar cod (Arctic Council, 2010; Stortini et al., 2016; Huserbråten et al., 2019) potentially leaving thermal niches open for other species to occupy (Laurel et al., 2016). Such changes could lead to dramatic changes to the fish assemblages accessible to subsistence fishers and commercial harvesters.

While physiological models are available for gadid species, they have yet to be integrated with climate projections to understand how warming ocean temperatures will influence suitable habitat for critical early life stages of marine fish in the northwest Atlantic. Here, we apply egg survival and juvenile growth models of three coastal species of gadids (Atlantic cod, Polar cod, and Greenland cod), with different thermal affinities and life histories, to ocean temperature projections. We illustrate how mechanistic physiological models can be used to predict habitat distribution shifts and compare vulnerabilities of species and life stages with changing ocean conditions. We also derived a habitat suitability index based on the combined surfaces of each metric to predict areas and climatology periods where thermal conditions were suitable for both life stages.

Materials and Methods

Study Area

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).

Spawning Habitat and Growth Potential Projections

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.

TABLE 1

Table 1. Cod development models for growth potential and egg survival.

Early Life Stage Habitat Suitability

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.

Distribution Shifts in Development Metrics

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:

C = \frac{\sum_{i = 1}^{N} w_{i} \cdot L a t_{i}}{\sum_{i = 1}^{N} w_{i}}

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).

Results

Marine waters off the east coast of Canada are likely to experience high levels of warming by the end of the century under the IPCC business as usual emissions scenario (RCP 8.5). In the north of the defined study area (≥55° N), from 2005 to 2100, median ocean temperatures are projected to increase by 1.9°C in the surface layer and 2.4°C in bottom layer during spawning season (December–March; Supplementary Figures 1A–D). In the southern portion of the study area, these increases will be greater with a projected increase of 4.3°C in the surface layer and 4.4°C in the bottom layer (Supplementary Figures 1E–H). During the first juvenile growth season (July–September), water temperatures are also projected to increase between 2005 and 2100. In the north of the study area, median ocean temperatures are projected to increase by 4.8°C in the surface layer and 2.4°C in bottom layer (Supplementary Figures 2A–D) while in the southern portion there are projected increases of 5.6°C in the surface layer and 4.1°C in the bottom layer (Supplementary Figures 2E–H). There was good agreement in the projected values for egg survival, growth potential and thermally suitable habitat between climate models (i.e., low signal to noise ratio). In general, model results had greater uncertainty in later periods (2051–2075 and 2076–2099), which is likely due to an increased spread in model predictions through time.

The response of cod egg survival to ocean warming will be species-dependent. Throughout the study area, median egg survival in Polar cod decreased from 86.2 to 65.7% survival by 2100 (Δ−20.5% median egg survival; Figures 1B, 2A) while median egg survival increased in Atlantic cod from 60.2 to 81.4% survival (Δ21.2% median egg survival; Figures 1B, 3A), and in Greenland cod from 34.3 to 53.3% survival (Δ19.0% median egg survival; Figures 1B, 4A). In all cases, these changes accompanied shifts in the distribution centroid and a northward shift in areas most suitable for egg survival (Figure 1A). For Polar cod, there was a poleward shift of 3.1° latitude in the distribution centroid by 2100 coupled with an ∼45% decrease in the area of spawning habitat (based on the 50^th percentile of baseline egg survival; Figure 5A). Remaining suitable spawning areas for Polar cod are predicted to be largely restricted to waters off the northern tip of Labrador (57°N–60°N) by the end of this century. Atlantic cod are forecasted to experience a smaller shift in their distribution centroid (1.2° north) by 2100, and egg survival is projected to increase by 53% in areal extent over this period. Atlantic cod spawning habitat is also forecasted to shift poleward from latitudes under 45°N to latitudes between 45°N and 50°N as warming at the southern edge of the study area exceeds the thermal optimum for egg survival. Lastly, like Polar cod, Greenland cod are forecasted to experience a large shift in their distribution centroid (3.1° north) by 2100. This is also due to decreased egg survival at the southern edge of the study area in conjunction with increased egg survival in the middle and northern portion of the study with warming. Similar to Atlantic cod, by 2100 there is a projected 27% increase in suitable spawning habitat area for Greenland cod, as poleward regions become more thermally favorable over current conditions.

FIGURE 1

Figure 1. Shift in the distribution centroid between the historical baseline period (1981–2005) and each future climatology period for egg survival (A), juvenile growth potential (C), and thermal habitat suitability index (E) for Polar cod, Atlantic cod, and Greenland cod. For each development metric, the distribution of individual raster cell values across the historic baseline period and future simulated periods are plotted (individual points; B,D,F) and summarized with boxplots showing the 25th, 50th, and 75th percentiles of plotted values for each species. Tails on boxplots represent 1.5 times the interquartile range.

FIGURE 2

Figure 2. Baseline and future projections of habitat suitability for Polar cod (B. saida) based on temperature-dependent egg survival (A), temperature-dependent juvenile growth potential (B), and geometric means of combined egg and juvenile habitat suitability (C). Dark gray cells indicate zero values while cells with high uncertainty (small change relative to spread across models) are indicated by a (+) symbol. The latitudinal component of the distribution centroid is denoted by a dashed gray line across each plot. Marginal plots along the y-axis show the mean cell values across latitude. The current distribution of B. saida is shown on the baseline habitat suitability (stippled area; bottom left panel).

FIGURE 3

Figure 3. Baseline and future projections of habitat suitability for Atlantic cod (G. morhua) based on temperature-dependent egg survival (A), temperature-dependent juvenile growth potential (B), and geometric means of combined egg and juvenile habitat suitability (C). Cells with high uncertainty (small change relative to spread across models) are indicated by a (+) symbol. The latitudinal component of the distribution centroid is denoted by a dashed gray line across each plot. Marginal plots along the y-axis show the mean cell values across latitude. The current distribution of G. morhua is shown on the baseline habitat suitability (stippled area; bottom left panel).

FIGURE 4

Figure 4. Baseline and future projections of habitat suitability for Greenland cod (G. macrocephalus) based on temperature-dependent egg survival (A), temperature-dependent juvenile growth potential (B), and geometric means of combined egg and juvenile habitat suitability (C). Gray cells indicate zero values while cells with high uncertainty (small change relative to spread across models) are indicated by a (+) symbol. The latitudinal component of the distribution centroid is denoted by a dashed gray line across each plot. Marginal plots along the y-axis show the mean cell values across latitude. The current distribution of G. macrocephalus is shown on the baseline habitat suitability (stippled area; bottom left panel).

FIGURE 5

Figure 5. Area (km²) by time period where cell values are greater than the 50^th percentile (threshold noted in bold after the last time point for each species) of baseline values (1981–2005) for egg survival (A), growth potential (B), and habitat suitability (C).

In contrast to egg survival, juvenile growth potential will likely decrease for all cod species by the end of the century due to warming. The median growth potential in the study area is projected to decrease from 1.0% mass⋅day^–1 to 0.4% mass⋅day^–1 (Δ−0.6% mass⋅day^–1 in median growth potential) in Polar cod (Figures 1D, 2B), from 3.6% mass⋅day^–1 to 2.4% mass⋅day^–1 (−Δ1.2% mass⋅day^–1 in median growth potential) in Atlantic cod (Figures 1D, 3B), and from 2.0% mass⋅day^–1 to 1.6% mass⋅day^–1 (−Δ 0.4% mass⋅day^–1 in median growth potential) in Greenland cod (Figures 1D, 4B). Juvenile Polar cod habitats are projected to shift 5.3° latitude to the north by 2100 (Figure 1C) as the areal extent of suitable juvenile growth habitats decline by 28% in the study region (Figure 5B). Further, by 2100 growth potential in most grid cells south of 50°N are predicted to fall under 0.5% mass⋅day^–1 (equivalent to ∼33% of maximum growth potential) for Polar cod. Atlantic cod are projected to see a large decrease (74%) in the total area of suitable growth potential habitat by 2100 alongside a shift in the distribution centroid of over 4.8° northward by 2100. Greenland cod are predicted to experience a slight decrease in median growth potential by the end of the century as the highest juvenile growth potential shifts from 45°N–50°N to areas north of 50°N. There is also a projected 63% decrease in the extent of suitable juvenile growth habitats for Greenland cod by 2100. The shift in distribution centroid for growth potential in Greenland cod largely mirrors the shift in the distribution centroid for Atlantic cod with only slight differences through time (shift 4.7° north by 2100).

An examination of HTS (i.e. combined egg survival and juvenile growth potential) suggests Polar cod will be negatively impacted by warming while HTS for Atlantic cod and Greenland cod will increase by the end of century. The median habitat suitability in the study region decreased from 0.70 to 0.37 by 2100 in Polar cod (Figures 1F, 2C) whereas habitat suitability in Atlantic cod remained relatively stable (from 0.53 to 0.51; Figures 1F, 3C) and Greenland cod habitat suitability increased from 0.28 to 0.40 (Figures 1F, 4C). The HTS metric was also forecasted to shift northward by 2100 in all three cod species (Figure 1E), with the largest shift occurring in Polar cod (4.9° latitude) followed by Greenland cod (4.8° latitude), and Atlantic cod (3.9° latitude). By the end of the century, there will likely be few areas of HTS habitat to support both eggs and juveniles below 50°N for Polar cod, as most grid cell values below that latitude approached zero. By 2100 there is a projected 35% decrease in the total area of habitat with HTS for Polar cod (Figure 5C). Initially there is an increase in the area of HTS habitat for Atlantic cod, but these gains were followed by sharp declines after 2050 to the end of century, with a projected 14% net decrease overall in areas of HTS. In contrast to the other gadid species, there is a projected 7% increase in areas of HTS in Greenland cod by 2100.

Discussion

Species-Specific Habitat Shifts

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.

Ecosystem Considerations

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.

Model Caveats and Assumptions

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.

Conclusion

The shifts in habitat distribution forecasted in this study will likely have far reaching consequences on ecosystem functioning and services (Pecl et al., 2017), which will in turn lead to socioeconomic impacts on communities and industries that rely on existing fisheries (Fulton, 2011; Hare et al., 2016; Ojea et al., 2020). The massive and extended collapse of the northern cod in the northwest Atlantic in the 1990s exemplifies the risks of not accounting for ecosystem change even for highly studied, regulated, and managed fisheries (Rose, 2005). Understanding the timing and magnitude of shifts in distribution and quality of habitat will help managers and resource users to factor projected trends into industry investment decisions and prepare for potential adjacency/historic access licensing conflicts. Although model projections from our study are unlikely to be helpful on an annual basis for management, they can inform longer term management decisions such as identifying acceptable levels of risk for harvesting and development (McHenry et al., 2019; Rose et al., 2000), MPA network design (McHenry et al., 2019; Brito-Morales et al., 2020) as well as monitoring, delineating stock boundaries, optimizing survey design, identifying vulnerable species, and projecting their critical habitat. Slowing the rate of carbon emissions is critical to minimizing ecosystem impacts and will clearly require action beyond fisheries science and management. Nevertheless, tools such as those implemented in this study can be used to understand the cost of inaction and to sustain the resilience of marine ecosystems in the interim.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

DC and BL conceived the study. CK and TH conducted the analysis with guidance from DC, BL, and TK. All authors contributed to writing the manuscript.

Funding

Funding was provided by the Nunavut Government and DFO’s Arctic Research Fund.

Conflict of Interest

TH is employed by Wood E&IS GmbH.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Thomas Hurst (NOAA) and two reviewers for their reviews of an earlier version of the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.764072/full#supplementary-material

References

Anderson, R. P. (2013). A framework for using niche models to estimate impacts of climate change on species distributions. Ann. N. Y. Acad. Sci. 1297, 8–28. doi: 10.1111/nyas.12264

PubMed Abstract | CrossRef Full Text | Google Scholar

Araujo, M. B., and Guisan, A. (2006). Five (or so) challenges for species distribution modelling. J. Biogeogr. 33, 1677–1688. doi: 10.1111/j.1365-2699.2006.01584.x