Suitability of foraging habitat for Eubalaena glacialis under future climate scenarios in the Northwest Atlantic Open Access

The North Atlantic right whale (Eubalaena glacialis) is one of the most critically endangered baleen whale species (listed as “Critically Endangered” on the red list of the International Union for Conservation of Nature (IUCN); Cooke, 2020), with fewer than 372 (+11/−12) individuals remaining in the population as of 2023 (Linden, 2024). E. glacialis, including the Northwest Atlantic population, was the target of intensive commercial whaling from the 12th century (Aguilar, 1986) until the practice was banned in 1935. The population was in slow recovery until 2010 (Waring et al., 2016), at which point the population once again decreased and sightings in traditionally used foraging habitats dropped (Davies and Brillant, 2019; Davies et al., 2019; Record et al., 2019; Pettis and Hamilton, 2024). Vessel collisions and fishing gear entanglements are significant sources of mortality for present-day populations: 85% of diagnosed deaths between 2010 and 2016 were attributed to entanglement and the other 15% to ship strikes (Kraus et al., 2016). Effective conservation strategies for E. glacialis often rely on knowledge of foraging habitat, such as lowering vessel speeds in historic foraging habitat (Laist et al., 2014). Shifting E. glacialis foraging habitat reduces the effectiveness of these strategies, and new, widespread mitigation measures are needed to ensure the survival of the species (Pettis and Hamilton, 2024).

Eubalaena glacialis inhabit the east coast of North America between Florida, United States, and Newfoundland, Canada (Winn et al., 1986; Kraus and Rolland, 2007). Calving occurs during the winter months near Florida and Georgia, and foraging has occurred historically in Cape Cod Bay during winter and in the Gulf of Maine and Bay of Fundy and on the Scotian Shelf in summer (Kraus and Rolland, 2007). Past management of E. glacialis had assumed implicitly that the species’ seasonal migration patterns would remain consistent (Record et al., 2019). However, many historical E. glacialis foraging habitats are experiencing some of the fastest rates of warming in the global ocean (Pershing et al., 2015). In 2010, E. glacialis populations shifted from their traditional foraging patterns, with a substantial proportion moving to occupy southern New England during the winter and the Gulf of Saint Lawrence in summer (Davis et al., 2017; Davies and Brillant, 2019; Simard et al., 2019; O’Brien et al., 2022). This change and the decrease in population size coincide with changes in the Gulf Stream and Labrador Current pathways that impacted bottom water temperatures and prey availability in the Gulf of Maine (Record et al., 2019; Friedland et al., 2020; Gonçalves Neto et al., 2021). Changes in E. glacialis distributions decrease the effectiveness of current conservation strategies relying on traditionally used foraging habitats, increasing the risk of anthropogenic mortality (Kraus et al., 2016; Record et al., 2019). As effective E. glacialis surveys are technically and logistically intensive, additional habitat identification methods will be useful to support survey area selection.

Eubalaena glacialis distribution during its main foraging season is traditionally linked to the abundance of its main prey species, Calanus finmarchicus, a lipid-rich copepod critical to the Northwest Atlantic ecosystem (Wishner et al., 1988; Murison and Gaskin, 1989; Mayo and Marx, 1990; Pershing and Stamieszkin, 2020). C. finmarchicus is an energy-rich prey resource for E. glacialis because later stages (C5s) accumulate high amounts of lipid to form dense layers during diapause (Michaud and Taggart, 2007; Ji, 2011). Historically, C. finmarchicus has occupied regions from the Irminger and Labrador seas near Greenland to the northern edge of the Mid-Atlantic Bight (Conover, 1988). The preferred depth and temperature requirements of diapause make C. finmarchicus especially vulnerable to climate change (Wilson et al., 2016), though local oceanography has helped the species persist in the Maine coastal current (Runge et al., 2014; Ji et al., 2017). Other Arctic calanoid copepod species such as C. hyperboreus and C. glacialis are also potentially suitable prey species for E. glacialis (Lehoux et al., 2020). These species currently inhabit waters from the Scotian Shelf northward (Head and Pepin, 2010) and may become increasingly important if E. glacialis distributions continue to move toward higher latitudes (Davis et al., 2017; Plourde et al., 2019; Davis et al., 2020; Meyer-Gutbrod et al., 2021; Meyer-Gutbrod et al., 2023). E. glacialis often appears to select foraging areas where Calanus spp. abundances are above a “critical threshold” satisfying bioenergetic requirements (Mayo and Marx, 1990; Kenney and Wishner, 1995; Plourde et al., 2019; Gavrilchuk et al., 2020). Examining current distributions of Calanus spp. and how they may shift with climate change will provide insight on changes to E. glacialis habitat.

A variety of species distribution models (SDMs) have been developed to examine present and future distributions of Calanus spp. (e.g., Ji, 2011; Reygondeau and Beaugrand, 2011; Chust et al., 2014; Grieve et al., 2017; Record et al., 2018; Lehoux et al., 2024). However, few have focused on Calanus spp. abundances relative to bioenergetic needs of Eubalaena glacialis (e.g., Pendleton et al., 2012; Ross et al., 2023). Additionally, many studies only consider the U.S. or Canadian regions of the contiguous northeast continental shelf (e.g., Grieve et al., 2017; Plourde et al., 2019; Ross et al., 2023; Plourde et al., 2024) or do not have sufficiently high spatial resolution to accurately model the complex oceanography around the Gulf of Maine (e.g., Reygondeau and Beaugrand, 2011; Chust et al., 2014). By considering multiple prey species and their availability to E. glacialis, merging abundance data from both U.S. and Canadian waters, and using a climate model downscaled for the region, we have aimed to provide a new perspective on changing E. glacialis foraging habitat over the entire Northwest Atlantic Shelf.

This study examines how dense Calanus spp. aggregations and the distribution of Eubalaena glacialis foraging habitats in the Northwest Atlantic may change in future decades under climate change. Additionally, we consider how E. glacialis suitable foraging habitat may shift relative to historical conditions and the influence of environmental covariates on suitability predictions. We built SDMs to predict E. glacialis foraging habitat suitability based on Calanus spp. abundance thresholds (Ross et al., 2023). Input prey aggregation data were obtained by comparing the historical abundance of dominant prey species C. finmarchicus and C. hyperboreus across the Mid-Atlantic Bight, Gulf of Maine, Scotian Shelf, Gulf of Saint Lawrence, and Newfoundland and Labrador shelves to a dry weight threshold. C. glacialis was omitted from this study because there was insufficient high-density data to allow reliable modeling. We then built boosted regression tree SDMs using the prey aggregation data and covariates obtained from a downscaled climate model for the region (Brickman et al., 2021). The SDM predictions were adjusted by basic E. glacialis bioenergetic correction factors (Plourde et al., 2019; Gavrilchuk et al., 2020) to obtain predictions of potential suitable E. glacialis foraging habitat. The models were used to project potential foraging habitats by coordinates, month, and prey species to 2046–2065 and 2066–2085 under the Intergovernmental Panel on Climate Change (IPCC) representative concentration pathway (RCP) 4.5 and 8.5 scenarios. The prey-specific projections were also combined to create comprehensive projections of potential new and lost E. glacialis foraging habitats over the entire Northwest Atlantic Shelf. The resulting habitat maps can be used to examine shifts in E. glacialis potential foraging habitat and inform new target areas for survey or conservation efforts.

The study domain extended from 35.2°N to 57.6°N and 77.0°W to 42.5°W (Figure 1), including the Mid-Atlantic Bight; Gulf of Maine; Scotian Shelf; Gulf of Saint Lawrence, further divided into the northern and southwestern regions; and Newfoundland and Labrador shelves. This domain includes the entire extent of the contributing Calanus spp. datasets as well as historically important Eubalaena glacialis foraging habitats, including Cape Cod Bay (Mayo et al., 2004) and the Great South Channel (Cetacean and Turtle Assessment Program, 1982; Kenney and Wishner, 1995) within the western Gulf of Maine, the Roseway Basin in the western Scotian Shelf, and the Bay of Fundy (Brown et al., 2014; Plourde et al., 2019). The domain also includes areas where E. glacialis have foraged since 2015, such as the Scotian Shelf and Gulf of Saint Lawrence (Plourde et al., 2019; Meyer-Gutbrod et al., 2023).

Environmental covariates were sourced as curvilinear point data from the Bedford North Atlantic Model (BNAM; Brickman et al., 2016). Climate trajectories were generated by BNAM by applying downscaling techniques to force ocean models of the North Atlantic with 1/12° spatial resolution (approximately 65 km², within the study area). Historical climate models were forced with the Coordinated Ocean–Ice Reference Experiments normal year (CNY) dataset and represent an average 1990–2015 climatology (Brickman et al., 2016; Wang et al., 2019). For future models for year-2055 (averaged 2046–2065 climatology) and year-2075 (averaged 2066–2085 climatology), RCP 4.5 and RCP 8.5 pathways were forced by adding IPCC circulation model anomalies to the CNY data. Monthly averaged covariates selected from the BNAM climate model for use were mixed layer depth, zonal and meridional current surface velocity (U and V, respectively), bathymetry (i.e., seabed depth), sea surface and sea bottom temperature (SST and SBT, respectively), and sea surface and sea bottom salinity (SSS and SBS, respectively). Calanus fitness has a complex interaction with the vertical structure of the water column and source waters (Plourde et al., 2024). We selected mixed layer depth, bathymetry, and surface/bottom temperature and salinity variables because of their key roles in these interactions. Velocity and mixed layer depth also have hypothesized roles in the formation of dense Calanus patches (Xue et al., 2022).

We obtained Calanus finmarchicus and C. hyperboreus integrated water column abundance data (ind m⁻²) from the National Oceanic and Atmospheric Administration (NOAA) Fisheries Ecosystem Monitoring Program (EcoMon) survey (US DOC/NOAA/NMFS) and the Fisheries and Oceans Canada (DFO) Atlantic Zone Monitoring Program (AZMP; Maillet et al., 2019; Casault et al., 2020; Blais et al., 2021; Figure 1). We chose to subset the available abundance data to the years 1990 through 2015 to correspond with the years represented by the modeled historical covariates. EcoMon surveys sampled zooplankton using a 333 µm mesh net to a maximum depth of 200 m; data have been collected since 1977 from the Gulf of Maine to the Mid-Atlantic Bight using stratified random sampling (Kane, 2007; Richardson et al., 2010). AZMP samples were obtained using a 202 µm mesh net through the entire water column; data were collected at set stations since 1999 on the Labrador, Newfoundland, and Scotian shelves and slope waters as well as in the Gulf of Saint Lawrence and Gulf of Maine (Therriault et al., 1998). EcoMon net tows that may underestimate zooplankton counts because of bathymetry deeper than 200 m were corrected using a generalized additive model (Text S1; Figure S1).

We considered only late-stage abundance data (i.e., C4–C6) for both Calanus finmarchicus and C. hyperboreus, where possible, because these lipid-rich stages are the primary target of Eubalaena glacialis foraging and were caught most effectively in the zooplankton sampling (Mayo et al., 2001; Plourde et al., 2019; Sorochan et al., 2019; Lehoux et al., 2020). C. hyperboreus was not a counted species in the EcoMon dataset. However, the zooplankton specimens identified as unidentified Calanus spp. in the EcoMon data were likely primarily C. hyperboreus and were therefore used for the C. hyperboreus model. The abundance counts for the unidentified Calanus spp. classification group were very low overall, aligning with expected low C. hyperboreus abundance in the area (Casault et al., 2020). While the unidentified Calanus spp. data are unstaged, there is evidence that C. hyperboreus reproduction primarily occurs further north in the Gulf of Saint Lawrence (Herman et al., 1991; Head et al., 1999). This suggests that the inclusion of potential sampled C1–C3 individuals in the Gulf of Maine would have little effect on abundance counts. For the C. hyperboreus dataset, we also removed 137 records from the AZMP dataset with high SBT (>5°C) along select transects and stations on the Scotian Shelf. These points likely represent outliers due to advective transport (Brennan et al., 2019).

We compared Calanus spp. abundance data to a critical water column biomass threshold, τ_b, to generate a binary presence/absence response variable for modeling. “Presence” was regarded as the presence of an aggregated abundance of large Calanus. In order to have a single threshold for both C. finmarchicus and C. hyperboreus while acknowledging differences in size, we converted raw abundance data to estimated water column biomass using species and region-specific estimates of individual dry weight from Sorochan et al. (2019; Table S1). Because the vertically corrected EcoMon datasets do not provide stage-specific abundances, we used C5 dry weight estimates to represent the average weight of late-stage individuals. C5 individual dry weight estimates were chosen as the most representative value because C5 size is between C4 and C6 sizes and, for C. finmarchicus, is the most lipid-rich and desirable stage for Eubalaena glacialis foraging (Sorochan et al., 2019). For each species and region, all water column abundance records (ind m⁻²) were multiplied by the species-specific individual dry weight to estimate water column biomass (g m⁻²). Off-shelf records had regions assigned by transect.

The critical biomass threshold, τ_b, is a static value based on the Calanus spp. water column biomass (g m⁻²) at which foraging should be energetically suitable for Eubalaena glacialis in waters shallower than 150 m in depth. Although some research has been conducted to observe and determine E. glacialis foraging thresholds (e.g., Mayo and Marx, 1990; Baumgartner and Mate, 2003; DeLorenzo Costa et al., 2006), most estimates are region-specific and vary due to differences in sampling and methods. A recent literature synthesis identified a range of 10,000–40,000 ind m⁻² as representative of reported E. glacialis foraging thresholds for C. finmarchicus in the Gulf of Maine (Ross et al., 2023). Random forest models using these abundance thresholds performed well at identifying E. glacialis foraging habitat and can reasonably translate to the larger area and comparable boosted regression tree algorithm of our study (Li and Wang, 2013). We chose to adopt an abundance threshold of τ_a = 30,000 ind m⁻², as it falls between the two thresholds determined by Ross et al. (2023) and maximizes ecological viability and model performance. The τ_a = 30,000 ind m⁻² value was multiplied by the individual dry weight estimate of 1.95 × 10⁻⁴ g ind⁻¹ for C. finmarchicus in the Gulf of Maine (Table S1) to convert the abundance threshold, τ_a, to a critical biomass threshold, τ_b, of 5.85 g m⁻² (summary of all terms in Table 1). This critical biomass threshold was therefore applied to our Calanus spp. water column biomass estimates.

To generate the covariates for the model, the Calanus spp. water column biomass data were matched spatially and temporally to the nearest BNAM historical climatological environmental covariates (Brickman et al., 2016). Points with Calanus spp. water column biomass above τ_b were considered a “τ_b-patch” and represent areas where Calanus spp. would be above minimum requirements for Eubalaena glacialis foraging (Table 1). The final C. finmarchicus and C. hyperboreus τ_b-patch datasets contained 19,222 and 21,344 samples, respectively. Positive τ_b-patch values made up 13.7% and 7.5% of the total points, respectively.

Cross-validated boosted regression tree models were trained on the Calanus finmarchicus and C. hyperboreus τ_b-patch datasets and environmental covariates. The models were designed for classification, meaning that they predicted presence/absence (τ_b-patch/no τ_b-patch) rather than a quantitative abundance of Calanus. We built separate models for each species to allow for species-specific responses to both the covariates (SDM training) and future environmental conditions (projections). Boosted regression trees were built using tidymodels version 1.1.0 (Kuhn and Wickham, 2020) and xgboost version 1.7.3.1 packages (Chen and Guestrin, 2016; Chen et al., 2024) within R 4.2.2 (R Core Team, 2023). Given a set of environmental covariates, the models predicted the probability of a τ_b-patch for their specified Calanus species on a scale from 0 to 1.

All available historical environmental covariates identified in Section 2.2 were used in the models, with the exception of SSS for the C. hyperboreus model. The SSS covariate was removed due to a strong but non-mechanistic association in the Gulf of Saint Lawrence that creates spurious results. U and V were transformed into a single current velocity magnitude covariate (velocity = sqrt[U² + V²]). Bathymetry was logarithmically base-10 transformed to reduce the impact of deep-water records. Additionally, all environmental covariates were normalized to have a mean of 0 and a standard deviation of 1. Boosted regression tree models are capable of handling non-normalized data (Elith et al., 2008), but normalization allowed for unbiased comparison of covariate importance using the vip version 0.3.2 package (Greenwell and Boehmke, 2020) in R 4.2.2. Covariate importance was calculated as the fractional contribution of each feature to model output. Month was included as a categorical explanatory variable to examine seasonality in modeled responses. We checked for spatial autocorrelation using the Moran’s I test in ape version 5.8.1 in R, in which values range from −1 to 1 and 0 represents a random distribution (Gittleman and Kot, 1990; Paradis and Schliep, 2019).

The boosted regression tree algorithm was chosen to maximize predictive performance and allow for highly detailed responses. This approach gives an advantage when predicting threshold-like behavior, allowing the algorithm to capture sharp jumps in the variable responses. In comparison with other algorithms, prediction onto independent data by boosted regression trees is also potentially less compromised by overfitting (Elith et al., 2008). Boosted regression trees for each species were tuned using tidymodels workflow sets by training and comparing models with varying hyperparameters to optimize the area under the receiver operating characteristic curve (AUC). AUC is a common method for evaluating SDMs and measures the model’s ability to distinguish between classes (Fielding and Bell, 1997). AUCs above 0.5 indicate better performance than a random model, and higher values are considered useful for estimating actual species distributions (Jiménez-Valverde, 2012). AUC was calculated using inbuilt tidymodels methods (Kuhn and Wickham, 2020). Models were highly insensitive to tuning; that is, model performance did not change significantly with different tuning settings (Elith et al., 2008). Hyperparameters were thus additionally chosen to minimize complexity and reduce overfitting. The tuned hyperparameters for both boosted regression trees set the number of trees at 500, tree depth at 4, learning rate at 0.1, random explanatory variables per tree at 5, and minimal node size at 10.

Cross-validation modeling techniques aggregate duplicate modeling algorithms trained on varying samples of the training dataset, improving overall performance and reducing the risk of overfitting (Berrar, 2019). In this study, we implemented the Calanus finmarchicus and C. hyperboreus boosted regression tree models by implementing Monte Carlo cross-validation with 100 subsamples. For each subsample, 75% of the dataset was used for training and 25% for testing, with the split stratified by τ_b-patch. The tuned boosted regression trees were trained on each of the 100 cross-validation data splits, resulting in 100 submodels for each species. The τ_b-patch probability predictions were calculated as the median value across all submodel predictions. The 95th percentile interquantile ranges for each prediction were calculated as the range between the 97.5th percentile and 2.5th percentile submodel predictions and serve as a measure of uncertainty.

After tuning and training, the True Skill Statistic was calculated using tidymodels (Kuhn and Wickham, 2020) as an additional method to verify model performance. The True Skill Statistic assesses the total proportion of true positives and ranges from −1 to 1, where scores above 0 correspond to better performance than a random model (Allouche et al., 2006). Monthly AUCs were extracted by calculating the median AUC across submodels for predictions in the testing datasets, grouped by month. The 95th percentile interquantile ranges for overall and monthly AUCs were calculated as the 2.5th and 97.5th percentile AUC of submodels. Finally, custom-built response curves were constructed to analyze model responses to variations in individual covariates (Text S2; Elith et al., 2005). The two trained, cross-validated, boosted regression tree models were then projected to the five climate scenarios to obtain 10 total projections of τ_b-patch probability: five per species, two per scenario.

The τ_b-patch predictions are based on a critical Calanus spp. water column biomass (g m⁻²) threshold used as a proxy of biomass density (g m⁻³) that would be suitable for Eubalaena glacialis foraging at depths shallower than 150 m. However, E. glacialis decreases daily ingestion time when prey is located deeper than 150 m (Baumgartner et al., 2017; see table in Gavrilchuk et al., 2020) and requires proportionally higher prey density to meet their energetic needs. Additionally, maximum Calanus spp. biomass density (g m⁻³) decreases with water column depth (Plourde et al., 2019). Therefore, the water column biomass (g m⁻²) required to maintain sufficient Calanus spp. density (g m⁻³) for E. glacialis increases with bottom depth. In order to account for these foraging considerations, we multiplied each projected τ_b-patch value by a “bioenergetic correction factor” calculated as a function of depth. We did not apply this correction factor earlier so that there would be a sufficiently high proportion of presences in the dataset for modeling. The bioenergetic correction factor function was obtained by multiplying functions of relative maximum prey density (from Plourde et al., 2019) and relative time spent foraging by a resting E. glacialis (derived from Gavrilchuk et al., 2020; Figure 2). We chose to consider resting E. glacialis foraging time (i.e., foraging time for an E. glacialis individual that is not pregnant or lactating) because it maximizes viable habitat in our results and is more representative of the “average” E. glacialis individual.

Each τ_b-patch prediction was multiplied by the bioenergetic correction factor for its associated bottom depth. This converted the τ_b-patch predictions into projections of suitable Eubalaena glacialis foraging habitats based on its basic foraging bioenergetic considerations (hereafter referred to as a “τ_h-patch,” where τ_h refers to the depth-varying critical threshold for potential E. glacialis foraging habitat; Table 1). “Habitat shift maps” for each future climate scenario, displaying new and lost potential suitable foraging habitats relative to historic distributions, were created by observing how τ_h-patch probabilities shifted between historical and future climate scenarios relative to a potential suitable habitat probability threshold, φ (Table 1). The threshold φ represents the minimum τ_h-patch probability required to consider an area potentially suitable for foraging. We explored and plotted multiple values of φ to consider the effect on projected shifts. We choose to focus on a φ value of 20%, as it falls roughly in the middle of the histogram distribution of τ_h-patch probability predictions in the testing dataset and maximizes ecological viability compared to other thresholds. Combined maps of τ_h-patch probability were taken by summing the probability of a τ_h-patch for either species at each 1/12° by 1/12° grid cell. Uncertainty for each prediction was calculated as the width of the 95th percentile interquantile range across submodels. We focus our conclusions on projections for the most extreme climate scenario, RCP 8.5 year-2075, representing no global efforts to mitigate carbon emissions, to highlight the strongest potential responses to climate change represented in the BNAM model.

Data were collected year-round for both surveys, albeit with varying effort by month (Figure 3a and b). The Calanus finmarchicus model had a median AUC of 0.881 (95th percentile interquantile range: 0.870–0.889). The C. hyperboreus model had a median AUC of 0.985 (95th percentile interquantile range: 0.981–0.988). The True Skill Statistic followed similar patterns to AUC and tended to be significantly above 0. AUC was consistently high month-by-month, with the exception of lower values in March–April for C. finmarchicus and July for C. hyperboreus (Figure 3c and d). The C. finmarchicus and C. hyperboreus modeled τ_b-patch probabilities correlated significantly to the underlying water column biomass data (Figure 3e and f; respective Spearman correlation coefficients of 0.75 and 0.64 with p < 2.2 × 10⁻¹⁶ for both). The C. finmarchicus model, while more correlated to the overall data, had significantly lower sensitivity and accuracy than the C. hyperboreus model (sensitivity [true positive rate]: 0.44 and 0.79, respectively; specificity [true negative rate]: 0.98 and 0.99, respectively; accuracy: 0.90 and 0.97, respectively). The Moran’s I test to assess spatial autocorrelation in the testing/training data detected statistically significant correlations of 0.108 and 0.563 for C. finmarchicus and C. hyperboreus, respectively (p < 0.05).

Relative covariate importance was similar between species (Figures 4 and S2). SBT was the most important variable for the Calanus hyperboreus model. Bathymetry was the most important variable for the C. finmarchicus model and second most important for C. hyperboreus. For both models, deep-water covariates were more important than sea surface covariates. Month and velocity contributed the least for both models.

The τ_h-patch response curves varied between species models (Figure 4). Response curves displayed minor high-frequency variations as well as larger low-frequency trends. Both models displayed a unimodal relationship to bathymetry, with the highest predictions for Calanus finmarchicus and C. hyperboreus observed at 170 m and 250 m, respectively. In contrast, the τ_b-patch response curves for bathymetry, reflecting model responses prior to the bioenergetic correction, both show a strong positive relationship. Both species had a strong negative relationship between mixed layer depth and τ_h-patch probability, with near-zero modeled probabilities at a mixed layer depth of 40 m and 30 m for C. finmarchicus and C. hyperboreus, respectively. The C. finmarchicus response curves depended unimodally on SST, SSS, and SBS, with respective peaks at 6°C, 25 ppt, and 34 ppt (Figure 4a). C. finmarchicus also displayed negative relationships with SBT and velocity, as well as a complex relationship to month with a significant peak during April. The C. hyperboreus response curves displayed unimodal relationships to SST and SBS, peaking at 0–4°C and 34 ppt, respectively (Figure 4b). The C. hyperboreus response curves also showed a threshold relationship to SBT, with near-zero probabilities predicted for temperatures above 5.5°C. The responses to velocity and month were mostly flat, with the exception of higher probabilities for near-zero velocities and summer months, respectively.

To describe model projections, we define predicted τ_h-patch probability as low where τ_h-patch probabilities were below 20%, moderate where probabilities were between 20% and 30%, and high where probabilities exceeded 30%. These cutoffs were chosen based on the distribution of probabilities and for consistency with φ.

In the combined Calanus spp. projections (“combined projections”), τ_h-patch probability at a point was calculated as the sum of τ_h-patch probability of C. finmarchicus and C. hyperboreus (Figure 5). C. finmarchicus τ_h-patch predictions predominantly explain observed areas of moderate suitability within the Gulf of Maine and on the western Scotian Shelf (Figure S3). C. hyperboreus τ_h-patch suitability areas were lower overall and concentrated in the Gulf of Saint Lawrence and on the Newfoundland and Labrador shelves (Figure S4). In regions with overlapped suitability, combined predictions were higher than either species predicted alone.

Model predictions varied monthly, with decreased τ_h-patch probability indicated between December and April (Figure 5). The combined projections suggested high τ_h-patch probability in the southwestern Gulf of Saint Lawrence in June and July. However, during winter months, τ_h-patch probability was generally higher in the Laurentian Channel—the central channel bisecting the northern and southwestern Gulf of Saint Lawrence—than the rest of the Gulf of Saint Lawrence. Projection output also suggested moderate τ_h-patch probability throughout summer in the Gulf of Maine, Mid-Atlantic Bight, Bay of Fundy, western Scotian Shelf, and Newfoundland and Labrador shelves. June and July, in general, displayed significantly higher τ_h-patch probabilities relative to other months, highlighting additional suitable habitat like the Flemish Cap region.

Projected uncertainty for combined predictions was generally proportional to the τ_h-patch prediction value, except for consistently low uncertainty for predictions in the Gulf of Maine (Figure S5). Uncertainty was highest in the southwestern Gulf of Saint Lawrence, particularly during June and July.

There were minor to moderate differences between each of the four future combinations of climate scenario and year (Figure S6). The greatest differences were observed between the RCP 4.5 year-2055 and RCP 8.5 year-2075 projections—the least and most extreme future projection scenarios, respectively. In general, later years and higher RCP scenarios decreased overall τ_h-patch probability. The most notable areas of τ_h-patch decreases from less to more extreme future projection scenarios were in the southern Gulf of Maine, southwestern Gulf of Saint Lawrence, and on the coastal Newfoundland and Labrador shelves. Some areas of moderate τ_h-patch probability increases between less and more extreme future projection scenarios were observed in the western Gulf of Maine and on the Scotian Shelf.

Habitat shift maps for RCP 8.5 year-2075 highlight potential areas of new, lost, and retained Eubalaena glacialis potential foraging habitat relative to a critical probability threshold φ of 20% (Figure 6). Winter months (November–April) showed little change (Figure S7). Overall, areas of lost potential habitat were greater than areas of gained potential habitat both across and within geographic regions (Figure S8). However, the relative proportion of retained habitat is dependent on φ (Figure S9). Lower φ values generally show more regions with suitable foraging habitat as well as a higher percentage of retained habitat, especially on the Newfoundland and Labrador shelves (Figure S9).

The combined projections showed a large loss of potential foraging habitat year-round in the central Gulf of Maine, particularly during the primary Eubalaena glacialis foraging months of July and August, although some potential habitat areas were gained along the western edge of the basin (Figure 6). Potential habitat was generally lost within the western Scotian Shelf, the southwestern Gulf of Saint Lawrence, and along the edges of the Laurentian Channel. Projections highlighted new potential habitat within the Roseway Basin in May and July, the Bay of Fundy in June, and in patches along the Newfoundland and Labrador shelves during June through August. Canadian waters, therefore, showed more new and retained potential habitat during summer months. Some areas of new and retained potential habitat were also projected along the edge of the Mid-Atlantic Bight in May and June.

Traditional conservation strategies for the critically endangered Eubalaena glacialis are challenged by recent distribution changes and would benefit from new, proactive approaches for management (Kraus et al., 2016; Davies and Brillant, 2019; Record et al., 2019). This modeling framework offers one approach for understanding future E. glacialis foraging habitat by combining and building on techniques from prior modeling studies: it focuses on a feeding-energetics-based threshold rather than direct prey abundance, it includes multiple prey species, and it projects future habitat scenarios across the full documented Canada–U.S. E. glacialis foraging domain. The model generally had high performance and projected potential changes in E. glacialis foraging habitat, which could inform survey effort and, indirectly, conservation recommendations.

Both the Calanus finmarchicus and C. hyperboreus models had strong performance based on AUC, accuracy, and sensitivity metrics (Pearce and Ferrier, 2000). C. finmarchicus AUC only dipped below 0.8 in March and April, and C. hyperboreus only dipped below 0.9 in April and July. These dips may be related to the spatial and temporal distribution in survey effort—for example, there are no AZMP surveys on the Scotian Shelf or in the Gulf of Saint Lawrence during July. The overall higher performance of the C. hyperboreus model is probably due to the strong and sharp threshold dependence on bottom temperature (Figure 4b). The high performance for both species during late summer and autumn correlates to key seasons for Eubalaena glacialis foraging, supporting the usefulness of the model (Murison and Gaskin, 1989; Pendleton et al., 2009; Brown et al., 2014). The low sensitivity of both models, particularly C. finmarchicus (0.44), suggests that the model may generate false negatives and thus our estimates may be overly conservative.

The most contributing covariates for both models were those correlated to deep-water conditions, that is, SBS and SBT. This finding is similar to previous modeling for Calanus finmarchicus (Grieve et al., 2017; Ross et al., 2023) and is consistent with the importance of deep-water diapause for both species. SBT influences diapause duration and, in turn, life cycle fitness (Varpe et al., 2009; Maps et al., 2014). SBS does not influence physiology significantly but is a strong indicator of water mass origin, which can determine copepod population supply (Ji et al., 2017; Record et al., 2019). Building Calanus SDMs based on surface oceanography alone is common (Pershing et al., 2009; Reygondeau and Beaugrand, 2011; Villarino et al., 2015), but the importance of deep-water variables in our model emphasizes the benefits of including subsurface oceanography.

Physics and water column structure influence the creation of Calanus patches. The strong contributions of bathymetry and, to some degree, mixed layer depth likely reflect how these covariates serve as proxies for oceanographic processes concentrating copepods. However, velocity had a weak contribution despite its more direct role in patch creation. Capturing the full physical process of patch formation likely requires finer-scale physics than our model, as well as a representation of copepod swimming against physics (Ji et al., 2012; Stegert et al., 2012). The strong contribution of bathymetry to the τ_b-patch model is likely mostly negated by the bathymetric correction, as the highest τ_b-patch predictions are from off-shelf locations, contrary to τ_h-patch.

The bioenergetic correction most affected the model response to bathymetry (Figure 4), as the highest τ_b-patch predictions and measured abundances are from off-shelf locations. Incorporating vertically resolved data (Krumhansl et al., 2018) would be one way to capture the seasonal variations in vertical distributions in models (Plourde et al., 2019; Plourde et al., 2024).

Thermal envelopes for high-density patches can differ from envelopes based on species abundance (Ross et al., 2023). The highest Calanus finmarchicus τ_h-patch probabilities were predicted at surface temperatures between 4°C and 8°C, which is similar to niches found in other modeling studies (Helaouët and Beaugrand, 2007; Reygondeau and Beaugrand, 2011). However, as in Ross et al. (2023), the bottom temperature range associated with high τ_h-patch probability was higher than in studies modeling abundance directly (Grieve et al., 2017). Reported thermal ranges for C. hyperboreus in the greater North Atlantic range from 0–2°C off Iceland to 5–7°C on the Norwegian shelf (Strand et al., 2020), falling within our predicted thermal envelope for surface temperature. The sharp C. hyperboreus response to bottom temperature was not seen in an earlier modeling study in this region (Albouy-Boyer et al., 2016); bringing out this feature probably required the broad Canada–U.S. dataset used here.

Overall, the model displayed generally high performance and displayed covariate responses consistent with prior studies and the Calanus life cycle.

Historical τ_h-patch projections for Calanus finmarchicus were in general agreement with historical Eubalaena glacialis sightings during summer foraging months. For example, the model predicted C. finmarchicus τ_h-patches during summer in the historically documented habitats of the Gulf of Maine and Bay of Fundy (Murison and Gaskin, 1989; Baumgartner and Mate, 2003). The training dataset does not contain measurements from the Bay of Fundy, so the success of the model in that habitat indicates some ability to interpolate in the near-shore zone. Additional highlighted areas on the Scotian Shelf and in the southwestern Gulf of Saint Lawrence have also been identified as suitable present-day foraging ground by previous modeling studies using C. finmarchicus prey layers (Plourde et al., 2019; Ross et al., 2023). The southwestern Gulf of Saint Lawrence, in particular, has had documented E. glacialis presence since 2010 (Simard et al., 2019; Meyer-Gutbrod et al., 2023). The high τ_h-patch probabilities south of New England and along the Mid-Atlantic continental shelf break are a reflection of high Calanus abundance measurements in the area and are consistent with increased Mid-Atlantic E. glacialis acoustic presence since the 2010 climate shifts (Davis et al., 2017; Murray et al., 2022). 

Historical projections for Calanus hyperboreus τ_h-patches, while outside of regions with historic Eubalaena glacialis survey effort, were generally associated with shallow areas known to have high C. hyperboreus abundance, such as the southwestern Gulf of Saint Lawrence and the Newfoundland and Labrador shelves (Head and Pepin, 2010). The highest C. hyperboreus τ_h-patch probabilities were predicted in the southwestern Gulf of Saint Lawrence, aligning with habitats predicted by Plourde et al. (2019). June had especially high predictions (65% positive τ_b-patch predictions during June, compared to 26% overall)—the absence of AZMP survey effort outside of the Gulf of Saint Lawrence in June may cause overestimation, although there is some evidence of denser, near-surface May and June aggregations in the lower Saint Lawrence Estuary and southwestern Gulf of Saint Lawrence (Plourde et al., 2003). Our model did not capture the possible high C. hyperboreus densities supplied by advection to the deeper regions of the Scotian Shelf (Sameoto and Herman, 1990; Herman et al., 1991), highlighting the uncertainty of a model not explicitly representing advection.

Projections of the two species together provide a basis for comparison of historical conditions against future projections. These projections should be regarded as scenario tests, not as forecasts. The τ_h-patch projections had only minor differences across future climate scenarios, which is consistent with model-based future predictions by Grieve et al. (2017) of Calanus finmarchicus abundance in that the main signal was an overall decrease in abundance for more extreme scenarios. However, because of the focus on threshold probability, our model contained finer-scale, varying distributions of increases and decreases.

We have focused here on the year-2075 RCP 8.5 projections for two reasons: the milder scenarios fall on a gradient between historical and year-2075 RCP 8.5 conditions (Figure S6), and observed rapid change in the Gulf of Maine region suggests that more extreme conditions could arrive earlier than expected (Pershing et al., 2015). In this scenario, nearly all significant areas of change are observed during the spring and summer, coinciding with the critical feeding months for Eubalaena glacialis. The large regions of loss projected within the Gulf of Maine have been projected by several other studies considering E. glacialis foraging habitats or Calanus finmarchicus abundance data (e.g., Reygondeau and Beaugrand, 2011; Grieve et al., 2017; Ross et al., 2021; Lehoux et al., 2024). These results provide further evidence that historical E. glacialis foraging habitat may continue to reduce under warming scenarios. Recent E. glacialis distribution shifts to areas such as the southwestern Gulf of Saint Lawrence may be driven more by the loss of historical prey aggregation areas than by the emergence of new ones (Meyer-Gutbrod et al., 2023). Our results, especially considering the disproportionate loss compared to gain of potential habitat, reinforce the idea that habitat loss could be a stronger driver of climate adaptation than habitat shift.

Overall, potential foraging habitat declined, but there were regions of new and retained habitat as well. In spring months (May–June), significant habitat was lost in the western Gulf of Maine, with a notable loss in the Gulf of Saint Lawrence in June, although with a small area of potential increase in the southwestern Gulf of Saint Lawrence in May. Potential areas of new habitat during these months included areas along various shelf breaks, as far ranging as the Flemish Cap and southern New England, as well as coastal areas around Newfoundland and Labrador and in the Bay of Fundy. Lost habitats in July–August were similar. The expanded areas on the western Scotian Shelf and in the Roseway Basin during this season align with predictions from Ross et al. (2021). Moving toward autumn (September–October) showed only small areas of potential habitat loss. The high habitat suitability areas that were retained and expanded along the Newfoundland and Labrador shelves and Flemish Cap are areas which have not generally been considered in prior Eubalaena glacialis modeling studies. The retained habitat on the Newfoundland and Labrador shelves, in particular, is consistent with the presumed stability in Calanus spp. abundance inferred by Plourde et al. (2024). The new and retained potential foraging habitats in the Bay of Fundy and Mid-Atlantic Bight also suggest that there still may be viable southern foraging for E. glacialis even with continued warming—consistent with some recent sightings. These results indicate areas where survey effort, when available, could be directed to give an early warning of right whale distributional shifts.

The methodology used in this study represented Calanus relative to an energetics-based feeding threshold, in contrast to directly modeling abundance (e.g., Plourde et al., 2019; Brennan et al., 2021; Plourde et al., 2024). By choosing to use a classification model, we focused on the information content most relevant to Eubalaena glacialis feeding at the expense of considering all possible abundances. Therefore, we cannot use the output of our model to generalize about Calanus spp. ecology, but rather about the associations between conditions and the probability of high-density Calanus measurements.

Both Calanus datasets contain sampling biases in space, time, and season. For example, while we cropped the data to match model covariates (1990–2015), AZMP sampling only began in 1999. We ran a test to examine the effect of this bias. Limiting both datasets to 1999–2015 left nearly all results unchanged except for lowered future τ_h-patch probability in July within the Gulf of Maine and on the western Scotian Shelf. Some areas within the study area are spatially undersampled, including historical Eubalaena glacialis foraging areas such as the Bay of Fundy and Cape Cod Bay (Murison and Gaskin, 1989; Mayo and Marx, 1990; Baumgartner and Mate, 2003). While the model displayed an ability to interpolate, some regions and months displaying extreme sampling bias are less reliable. Some degree of spatial autocorrelation was also present for both species. C. hyperboreus displayed higher autocorrelation, likely due to the strong latitudinal gradient for the species. Spatial autocorrelation, while a common source of uncertainty in SDMs, may lead to underestimated uncertainty and biased parameter estimates (Guélat and Kéry, 2018). The strong influence of variables with latitudinal variation, like SBT in the C. hyperboreus model, is likely affected by spatial autocorrelation and should be interpreted with caution.

As with any model, some ecological parameters in the model were approximated with constant or averaged values. Calanus dry weights (Plourde et al., 2019; Sorochan et al., 2019; Helenius et al., 2024) and stage composition, as well as Eubalaena glacialis foraging thresholds and target species (Plourde et al., 2019; Meyer-Gutbrod et al., 2023) depend on changing environmental conditions, regional variation, and internal physiological states. The sensitivity of E. glacialis SDMs to the τ_b foraging threshold was explored by Ross et al. (2023) and appears to primarily affect the scale of predictions rather than relative probabilities between regions. However, Ross et al. (2023) focused on foraging thresholds within the Gulf of Maine without considering variation in other regions. Lehoux et al. (2024) explored variation in Calanus dry weight and predicted that dry weight would decrease across the board by the end of the century, particularly within the Gulf of Saint Lawrence. This finding suggests that our approach using a constant dry weight may overestimate suitability. Our model could be improved by approximating these parameters with dynamic, region-based functions. The habitat shift maps are strongly dependent on our choice of φ (Figure S9), although general areas of loss and gained suitability are retained. Further research on the minimum viable φ for E. glacialis foraging would improve the reliability of retained suitability predictions.

Despite efforts to reduce overfitting via cross-validation, model selection, and hyperparameter choices, the high-frequency variations in the response curves (Figure 4) still suggest a degree of overfitting. Such overfitting is indicative of the model compensating for “data noise” around lower-frequency variations and likely does not correlate to mechanistic relationships (Merow et al., 2014). Because the amplitude of the high-frequency variations is much lower than that of low-frequency variations, overfitting is unlikely to have a major impact on our projections. However, it may result in overestimation of performance metrics and artificial small-scale variation in predictions.

Additional variables omitted from this model may also contain information important for modeling. Other non-Calanus copepod prey taxa, such as Centropages and Pseudocalanus (Mayo and Marx, 1990; Pendleton et al., 2009; Hudak et al., 2023), may also influence Eubalaena glacialis climate adaptation by contributing to additional suitable foraging areas. Our model does not consider potential Calanus adaptations to change (Kvile et al., 2018), although they may be limited (Hinder et al., 2014); these adaptations are better captured in fully coupled or individual-based models. Additional variables such as ocean pH may also influence Calanus spp. patches (Cripps et al., 2014; Thor et al., 2016; Fields et al., 2023). Finally, our model does not capture seasonal variation in vertical distribution (Plourde et al., 2019; 2024), which could be addressed by further research incorporating vertically resolved data (Krumhansl et al., 2018). These dependencies are currently common sources of uncertainty in Calanus models and are under ongoing research.

Environmental covariates were sourced from a downscaled physics model, which itself approximates reality and is subject to uncertainty. BNAM is a highly conservative climate projection for this region, and the area has already come close to the projected levels of warming from the model (Pershing et al., 2021). Given the disproportionate rate of warming in this region (Pershing et al., 2015), expecting these projected new habitats to emerge in a much nearer time frame than 2055 or 2075 is reasonable.

Finally, a prediction of suitable foraging habitat does not necessarily imply whether whales will occupy that space. The computed shifts in foraging habitat (Figure 6) illustrate the possible climate deficit experienced by whales, but there is also an adaptation deficit where whales must adapt to new conditions (Pershing and Pendleton, 2021). If the environment is in perpetual change, when and if Eubalaena glacialis will locate suitable foraging habitat before it shifts may be hard to determine. Further research on how whales make foraging decisions may help address this knowledge gap. Regardless, these projections offer an approach for targeting survey efforts toward areas most likely to support E. glacialis foraging.

Making projections into future climate scenarios has become a common practice in contemporary ecology, but how to best interpret and use these predictions is not always clear. Many decades may pass before predictions can be evaluated against data in terms of their benefits or harms to decision making. There is also the potential for rapid unexpected changes to render the predictions meaningless (Record et al., 2023), such as an early abrupt collapse of the Gulf Stream (Ditlevsen and Ditlevsen, 2023). The Calanus patch maps provided here are intended to provide one tool among many for responding to the changing ocean. Policy recommendations should not be made on models alone, but models can be used to help guide cost-limited survey efforts as we try to anticipate new foraging areas for right whales. Survey efforts, in turn, can inform conservation actions like establishing protected areas, seasonal regulations, and adaptive management. Given that some of the conditions predicted in the BNAM model are already being observed, we argue that these Calanus threshold projections can be put to use already and assessed against new whale-feeding locations as they are identified. We also recognize that this problem is complex and dynamic and these projections should be used alongside modeling efforts representing other important dynamics (e.g., Ross et al., 2021; Plourde et al., 2024), as well as with expert human judgment, to more fully inform decisions.

EcoMon data is accessible through the National Oceanic and Atmospheric Administration website (Northeast Fisheries Science Center, 2019). AZMP data as well as the oceanographic covariates from the BNAM climate model are publicly available through the Department of Fisheries and Oceans Canada (https://www.dfo-mpo.gc.ca/science/data-donnees/azmp-pmza/index-eng.html). Data requests should be coordinated through these agencies. Modeling code can be retrieved at https://github.com/BigelowLab/predicting_NARW.

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

Text S1–S2. Figures S1–S9. Table S1. Docx

Thanks to many people involved in collection and maintenance of the EcoMon and AZMP datasets. Additional thanks to David Brickman for developing the BNAM climate model.

Funding was provided by NSF grants #1849227 and #2307754, NOAA IOOS grant #NA24NOSX012C0020, and Bigelow Laboratory institutional funds.

There are no competing interests.

Contributed to conception and design: OHJ, SP, CL, CHR, CDO, NRR.

Contributed to acquisition of data: SP, CDO, HJW.

Contributed to analysis and interpretation of data: OHJ, SP, CL, BT, NRR.

Drafted and/or revised the article: OHJ, SP, CL, CHR, BT, CDO, HJW, NRR.

Approved the submitted version for publication: OHJ, SP, CL, CHR, BT, CDO, HJW, NRR.

How to cite this article: Johnson, OH, Plourde, S, Lehoux, C, Ross, CH, Tupper, B, Orphanides, CD, Walsh, HJ, Record, NR. 2025. Suitability of foraging habitat for Eubalaena glacialis under future climate scenarios in the Northwest Atlantic. Elementa: Science of the Anthropocene 13(1). DOI: https://doi.org/10.1525/elementa.2024.00044

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

Knowledge Domain: Ocean Science