Climate change is altering the marine environment at a global scale and these changes could affect the distribution and migration patterns of marine species throughout their annual cycle. Arctic regions are already experiencing some of the most dramatic changes in marine climate and there is a need for predictive models to understand how these changes could alter the spatio-temporal distributions of Arctic marine species. We used a species distribution model to predict potential future changes in the non-breeding distribution of thick-billed murres (Uria lomvia) from a colony in Hudson Bay, Canada, from 2021 to 2100 using three Coupled Model Intercomparison Project Phase 6 (CMIP6) climate scenarios: low (SSP1-2.6), intermediate (SSP2-4.5), and high (SSP5-8.5) emissions. Under the intermediate and high emissions scenarios, suitable habitat within Hudson Bay would become available year-round during the next century. This could lead to a portion of this migratory population becoming year-round residents within the next 80 years. We predicted a significant northward shift in the winter range, such that little or no habitat would be available below 55ºN by 2100. This shift would have significant implications for the murre harvest in Canada because the winter distribution would no longer include coastal Newfoundland where most harvesting occurs, particularly if murres from other colonies show a similar shift in distribution. Although there were projected changes in seasonal distributions under all three climate scenarios, dramatic re-distribution of non-breeding habitat could be avoided with policies that limit future emissions.

We used a previously published species distribution model to predict non-breeding distributions of thick-billed murres under different climate change scenarios (Patterson et al. 2021, Figure 1). The model was developed using geolocator tracks from 91 adult thick-billed murres tracked from a breeding colony on Coats Island, Nunavut (62.95ºN, 82.01ºW), collected over four annual cycles (2007-09, 2017-2019) (Patterson et al. 2021b). Full details on geolocator unit deployment and location estimates are outlined in Patterson et al. (2021). Location estimates were derived from light-level data using the R packages ‘TwGeos’ (Lisovski et al. 2016) and ‘probGLS’ (Merkel et al. 2016). Pseudo-absences representing habitat available to murres were randomly generated for ocean areas within a 200-1000 km envelope from all observed locations within the month of tracking at a 1:1 ratio of pseudo-absences to observations. Environmental conditions at observed locations were compared to environmental conditions at pseudo-absence locations paired by date. The model used daily time-varying environmental variables for sea surface temperature (SST, NOAA high-resolution SST, NOAA/OAR/ESRL PSL, https://psl.noaa.gov/), air temperature (NOAA Physical Sciences Laboratory NCEP/NCAR Reanalysis 1 https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.pressure.html), sea ice concentration (NOAA high-resolution ice cover, NOAA/OAR/ESRL PSL, https://psl.noaa.gov/), and wind speed (NOAA Physical Sciences Laboratory NCEP/NCAR Reanalysis 1 https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.pressure.html), as well as fixed predictors for bathymetry (ETOPO1 Global Relief Model, www.ngdc.noaa.gov/mgg/global/), slope, distance from colony, and day of year. All predictor variables were standardized to a common 0.25º x 0.25º raster resolution. The SDM was fit with a random forest, using the ‘ranger’ package (Wright & Ziegler 2017), to estimate probability of occurrence as a function of environmental conditions throughout the non-breeding season (Sept 1-May 31). Distance from colony and SST had the greatest influence on murre distributions (see Fig 2, Patterson et al 2021). From this fitted model we predicted expected distributions under different future environmental conditions (see below). A complete description of original SDM development is provided in Patterson et al. (2021).

All models participating in CMIP6 run a historical experiment using observed values for greenhouse gasses and other climate drivers over the period 1850-2014 (Eyring et al. 2016). These historical simulations are used to assess how well models simulate climate, document model characteristics, and ensure continuity across phases of CMIP (Eyring et al. 2016, O’Neill et al. 2016). We used the historical simulations from nine candidate models to determine which models gave credible historical distributions relative to our SDM predictions from the period for which remote sensing observations are available (1982-2014). We considered these nine CMIP6 models: ACCESS-CM2 (Dix et al. 2019), CanESM5 (Swart et al. 2019), CMCC (Lovato & Peano 2020a), EC-Earth3 (EC-Earth Consortium 2019), MIROC6 (Tatebe & Watanabe 2018), MPI-ESM (Wieners et al. 2019d), MRI-ESM2-0 (Yukimoto et al. 2019a), NorESM2-LM (Seland et al. 2019), and NorESM2-MM (Bentsen et al. 2019). These models were chosen because each predicted four relevant environmental variables (SST, air temperature, sea ice concentration, windspeed) at a daily frequency at a nominal spatial resolution of 100 km. CMIP6 data were obtained from https://esgf-node.llnl.gov/search/cmip6/.

We used the SDM model described in Section 2.1 to make predictions from 1982-2014 using the historical simulations from each model as inputs for SST, air temperature, sea ice concentration, and windspeed. SDM predictions were made at 7-day intervals for the non-breeding portion of the annual cycle (Sep 1 - May 31). Within each stage of the annual cycle (moult: 245-307 day-of-year (DOY); fall migration: 308-362 DOY; winter: 363-88 DOY; spring migration: 89-152 DOY) we calculated the median probability of use for each cell with predictions made from remote sensing observations and predictions made from the nine climate models. We used Schoener’s D (Schoener 1968, Warren et al. 2008) to compare the similarity between predictions based on remotely sensed data and each climate model. Schoener’s D was calculated as: D = 1 - 0.5 * ∑(|P₁ - P₂|) where, P₁ is predicted probability of use values from remotely sensed historical data and P₂ is predicted probability of use from the historical experiments of climate models. This statistic can take values between 0 (no overlap) and 1 (perfect overlap). Only raster cells with a value ≥ 0.5 for P₁ or P₂ were included in comparisons calculations. Cells with low predicted values in both the observed and predicted data were excluded to avoid having areas with a low probability of use strongly influencing the assessment, however model rankings were not sensitive to this threshold (range tested: 0.3-0.7). For each CMIP6 model, we calculated the mean Schoener’s D across the four annual cycle stages. Models with mean Schoener’s D greater than 0.90 were considered to adequately predict historical distributions. Predictions from models that met this criterion were averaged to create a composite model. Schoener’s D from the composite model was compared to individual models to determine whether the composite model provided better consistency with predictions from remote sensing data.

The CMIP6 provides predictions of environmental change under a range of predicted future climate scenarios, known as the Shared Socio-economic Pathways (SSP) (O’Neill et al. 2016). SSP1-2.6 encompasses the low end of future emissions considered within CMIP6 scenarios (O’Neill et al. 2016); this scenario assumes low challenges to climate mitigation and adaptation (O’Neill et al. 2017, Riahi et al. 2017). SSP2-4.5 is the ‘Middle of the Road’ pathway assuming intermediate challenges to mitigation and adaptation, consistent with historical trends observed over the past century (O’Neill et al. 2017, Riahi et al. 2017). Finally, SSP5-8.5 encompasses the high end of future emissions (O’Neill et al. 2016). This scenario of high future emissions is only feasible under SSP5, the ‘Fossil-fueled Development’ pathway, which assumes rapid growth of the global economy coupled with intensive fossil fuel development, high challenges to climate mitigation, and low challenges to climate adaptation (O’Neill et al. 2017, Riahi et al. 2017). Data from the SSP 1-2.6, SSP 2-4.5, and SSP 5-8.5 experiments were used to predict potential future distributions using the climate model with the highest Schoener’s D. Predictions were made at 7-day intervals for the non-breeding portion of the annual cycle (Sep 1 to May 31) for every year from 2020-2100.