Northwest range shifts and shorter wintering period of an Arctic seabird in response to four decades of changing ocean climate

: Climate change is altering the marine environment at a global scale, with some of the most dramatic changes occurring in Arctic regions. These changes may affect the distribution and migration patterns of marine species throughout the annual cycle. Species distribution models have provided detailed understanding of the responses of terrestrial species to climate changes, often based on observational data; biologging offers the opportunity to extend those models to migratory marine species that occur in marine environments where direct observation is difficult. We used species distribution modelling and tracking data to model past changes in the non-breeding distribution of thick-billed murres Uria lomvia from a colony in Hudson Bay, Canada, between 1982 and 2019. The predicted distribution of murres shifted during fall and winter. The largest shifts have occurred for fall migration, with range shifts of 211 km west and 50 km north per decade, compared with a 29 km shift west per decade in winter. Regions of range expan-sions had larger declines in sea ice cover, smaller increases in sea surface temperature, and larger increases in air temperature than regions where the range was stable or declining. Murres migrate in and out of Hudson Bay as ice forms each fall and melts each spring. Habitat in Hudson Bay has become available later into the fall and earlier in the spring, such that habitat in Hudson Bay was available for 21 d longer in 2019 than in 1982. Clearly, marine climate is altering the distribution and annual cycle of migratory marine species that occur in areas with seasonal ice cover.


INTRODUCTION
Climate change is altering the marine environment worldwide (Hoegh-Guldberg & Bruno 2010), changing the phenology and distribution of marine flora and fauna . Arctic surface water temperature increased at a rate of 0.5°C decade −1 from 1982 to 2017 (Meredith et al. 2019), and surface air temperature has increased twice as fast as the global average in the last 2 decades (Meredith et al. 2019), leading to the rapid loss of sea ice, an important physical stratum for wildlife. Consequently, many Arctic marine species are experiencing changes in distribution, abundance, and pheno logy, either as a direct response to physical changes in their habitat or indirectly through trophic interactions (Sydeman et al. 2015). For highly mobile, pelagic species, our knowledge of their reaction to habitat changes has been limited by our ability to observe animals at sea, especially outside of the breeding season. Given the strong seasonality in the Arctic, there is an urgent need to measure how habi-tat use and phenology change through the annual cycle, to understand how climate change is affecting Arctic marine life.
Ice cover and ocean temperature are important factors determining the large-scale distributions and abundance of marine species (Perry et al. 2005, Post et al. 2013. Ice directly affects polar marine mammals and marine birds by either facilitating or restricting access to prey (Tynan et al. 2009). Changing sea surface temperature (SST) and ice cover in the northwest Atlantic are associated with changes in the growth of Atlantic salmon Salmo salar (Friedland & Todd 2012), anadromous Arctic charr Salvelinus alpinus (Michaud et al. 2010), capelin Mallotus villosus (Carscadden et al. 2001), and Atlantic cod Gadus morhua (Drinkwater 2005), and, farther north, Arctic cod Boreogadus saida and zooplankton blooms (Welch et al. 1992, Beaugrand et al. 2003, Darnis et al. 2012. Moreover, seasonal sea-ice dynamics also play an important role in the timing of spring phytoplankton blooms, which can be a key factor at the end of the non-breeding season (Coppack & Both 2002, Søreide et al. 2010, Gaston et al. 2011, Leu et al. 2011, Post et al. 2013. Ocean warming could also have significant effects on the survival and reproductive success of many polar seabird species, through changes in the distribution, abundance, and availability or their prey (Croxall 2002, Sydeman et al. 2015.
Thick-billed murres Uria lomvia (hereafter, murres) are abundant and widespread Arctic seabirds with a circumpolar distribution. The species is considered an important indicator of Arctic marine ecosystems (Mallory et al. 2006, Barry et al. 2010, Michel et al. 2012. Through much of their range, murres migrate away from breeding areas as ice forms in winter, and return as ice recedes in spring. A longer ice-free period could influence the timing and extent of migration by murres, allowing them to remain within their breeding range longer. The decline of sea ice could also affect the availability of ice-associated prey (Hop & Gjøsaeter 2013), or the timing of peak prey availability relative to key periods of the murre annual cycle, such as chick-rearing. Increasing ocean temperature could affect the distribution, abundance, and size of their prey (Carscadden et al. 2001, Drinkwater 2005, von Biela et al. 2019, while simultaneously increasing competition with predatory fish species (Holsman & Aydin 2015). This type of complex trophic interaction has been suggested as the cause of mass mortality and breeding failure for common murres U. aalge in the northeast Pacific in response to an extended marine heatwave (Piatt et al. 2020). Increased frequency and intensity of storms could increase foraging costs for murres. Many seabird species have been shown to spend more time foraging or have lower feeding rates during inclement weather (Finney et al. 1999, Daunt et al. 2006, Elliott et al. 2014, and winter mortality events have been associated with periods of high wind (Harris & Wanless 1996, Frederiksen et al. 2008.
Species distributions models (SDMs) are important tools for predicting the current, past, and future distributions of wildlife (Elith & Leathwick 2009, Dambach & Rödder 2011, Robinson et al. 2011, Guisan et al. 2013 and have been used widely in terrestrial ecology for the last 30 yr. They have been used less in marine ecology (Robinson et al. 2011), most commonly for fish and marine mammals (Dambach & Rödder 2011). The proliferation of tracking studies on seabirds provides an opportunity to use SDMs to expand our understanding of habitat use of marine birds at sea . Species with high dispersal ability, such as seabirds, are more capable of tracking climate changes than more sedentary species, facilitating modelling of shifts in distribution (Araújo & Pearson 2005). Data from tagging studies provide continuous sampling of a species' distribution and habitat preference with less spatial and temporal bias than other visual surveys, especially in remote habitats where direct observation is difficult (Dambach & Rödder 2011. We used SDM and global location sensor (GLS) tracking data to model habitat use and examine past changes in the non-breeding distribution of murres from Coats Island, Nunavut, Canada. An SDM was developed for the non-breeding period (September to May) using tracking data collected over 4 nonbreeding periods (2007/08, 2008/09, 2017/18, and 2018/19). We used climate and physical oceanography variables to model the non-breeding distribution of murres. From this model, we predicted the historical distribution of murres from 1982 to 2019 using remotely sensed climate data. We used these predictions to map murre distributions during 4 non-breeding stages of the murre annual cycle (moult, fall migration, winter, and spring migration) and to test for long-term changes in these distributions. We expected to find that non-breeding distributions have shifted north as a result of warming ocean temperatures and declining sea ice cover, which are known to be occurring within the range of this population. We also tested for changes in the phenology of habitat availability within Hudson Bay. We expected that more habitat would be available for murres within Hudson Bay in fall and spring, due to declining seaice cover in Hudson Bay. Exploring the extent, mag-nitude, and direction of these past changes is an important first step to understanding how sensitive this species will be to future climate change.

GLS tracking
Tracking was conducted at Coats Island, in Hudson Bay, Nunavut, Canada (62.95°N, 82.01°W), a colony of 30 000 breeding pairs of thick-billed murres . As part of earlier tracking studies (Gaston et al. 2011, McFarlane Tranquilla et al. 2013, 3 types of geolocators (British Antarctic Survey), namely Mk5 (3.6 g), Mk7 (3.6 g), and Mk13 (1.8 g), were deployed at Coats Island in 2007 (n = 20) and 2008 (n = 20). All Mk loggers recorded maximum light levels at 10 min intervals, and a subset of these loggers also recorded temperature at 10 min intervals. In 2017 (n = 48) and 2018 (n = 45), we deployed LAT2800S geolocator− temperature− depth recorders (Lotek; 36 mm × 11 mm × 7.2 mm, 5.5 g) at the same colony. LAT 2800 loggers were programmed to collect light level, temperature, depth, and wet/dry state at 10 s intervals. All loggers were deployed during the summer on breeding adults captured on the nest, using a noose pole, while attending an egg or chick. Loggers were retrieved and data downloaded 1 to 2 yr later, during subsequent breeding seasons.

Location estimates
For LAT2800S loggers, we summarized maximum-recorded light levels at 5 min intervals prior to estimating twilight. Twilight was estimated using the threshold method in the 'TwGeos' package (Lisovski et al. 2016). We defined 2 behavioural modes: flying and on water. Flying was defined as any period where the sensor was dry and tag temperature was < 5°C; this temperature threshold was used to prevent periods of leg-tucking from being falsely classified as flying (Linnebjerg et al. 2014). Location estimates were calculated using a probabilistic algorithm with the 'probGLS' package in R (version 0.9.5, Merkel et al. 2016). The 'probGLS' method estimated 2 locations daily, at sunrise and sunset. At each time step, 1000 random particles were generated within the defined study area based on the observed twilight, random solar angles between −6° and −1°, and twilight error following a log-normal distribution (shape = 2.49, scale = 0.94). Using 'probGLS', we also incorporated additional information about habitat use and murre behaviour by weighting each random particle based on a land mask, sea ice cover, SST, and movement speed. Because murres do not use, or travel over, land during the non-breeding period, random particles over land received a weight of 0. Because murres cannot re main in areas with complete ice cover, random particles with greater than 90% ice cover (NOAA high-resolution ice cover, NOAA/OAR/ESRL PSL, https:// psl.noaa.gov/) were also assigned a weight of 0. For loggers with a temperature sensor (LAT 2800 and MK5), particles were weighted according to the similarity between re motely sensed SST and internal logger temperature (NOAA high-resolution SST, NOAA/ OAR/ESRL PSL, https://psl.noaa.gov/; Reynolds et al. 2007). Finally, random particles were weighted ac cording to a movement model limiting the distance travelled be tween consecutive locations based on realistic movement rates for murres; the movement model used different parameters for loggers with wet/dry sensors (LAT2800) that could estimate time spent in flight (Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m6790 p163 _ supp .pdf). One particle was randomly selected from the possible particles based on the assigned weights. These steps were repeated at each time step until an entire track was generated. The process was repeated to generate 100 possible tracks for each deployment. The most probable track was calculated as the geographic median of possible locations at each time step, and this track was used in mapping and estimates of migration timing. Full details of the parameters used in the probabilistic algorithm are provided in Table S1. We present maps of estimated tracks for each year, based on the most probable tracks estimated above (see Fig. 1).
To compare the timing of migration across tracking years, we calculated the latest date when each bird crossed 70°W in fall and spring. Murres from Coats Island migrate through Hudson Strait, and 70°W represents the halfway point of movement through this corridor. We used mixed effects models to test for differences in migration timing among years, with individual identity as a random factor to account for murres tracked over 2 yr. Mixed effects models were constructed using the 'lme4' package, version 1. 1-27 (Bates et al. 2015). Residual and q-q plots were used to check assumptions. The minimum and maximum dates when 95% of tracked birds migrated across years were used to summarize the fall and spring migration stages, respectively.

SDMs
We developed an SDM for the non-breeding period (September to May). SDMs assume that species are at an equilibrium with their environment and that all relevant environmental gradients have been sampled (Elith & Leathwick 2009). Tracking of smaller species is often limited to using GLS devices, which have lower spatial accuracy than other tracking methods (Phillips et al. 2004). Using paired deployments of satellite platform terminal transmitters (PTTs) and GLS loggers on black-browed albatrosses Thalassarche melanophris, Quillfeldt et al. (2017) found that device type (PTT versus GLS) had less influence on SDM accuracy and overlap in predicted distributions than the choice of SDM algorithm.
We considered 7 predictor variables in the model: bathymetry, slope, distance from colony (distance), day of year (DOY), SST, air temperature (air), sea ice cover (ice), and wind speed (wind). The 3 static environmental variables, i.e. bathymetry, slope, and distance from colony, were included in the model because they are likely biologically relevant to the species' distribution (Stanton et al. 2012). Bathymetry and slope are both likely to influence the distribution of prey, even with changing marine climate. Murres are primarily constrained to areas close to the colony during moult and spring migration; therefore, distance from the colony was included to ensure this pattern was included in the model. DOY was in cluded as a temporal variable to allow habitat preferences to change through the non-breeding period. Other static variables, i.e. longitude, latitude, and day length, were not included because we were interested in how the species responds to variation in climate (Stanton et al. 2012). Bathymetry data used the ETOPO1 Global Relief Model (www. ngdc. noaa.gov/ mgg/ global/). Slope was calculated from the bathymetry layer using the terrain function in the 'raster' package in R (Hijmans & Van Etten 2016). Daily-mean SST and ice were obtained from the European Space Agency Reprocessed Sea Surface Temperature Analysis (Merchant et al. 2019). Daily mean air temperature (2 m) and surface wind speed were obtained from NOAA Physical Sciences Laboratory NCEP/ NCAR Reanalysis 1 (https://psl. noaa. gov/data/ gridded/ data.ncep.reanalysis. pressure .html). Bilinear interpolation was used to resample environmental variables to a standard 0.25°s patial resolution, using the 'raster' package (Hijmans & Van Etten 2016).
Our SDM compared environmental predictors at observed locations from murre GLS tracks to available environmental conditions at pseudo-absence locations (Barbet-Massin et al. 2012). To incorporate uncertainty in GLS location estimates into the modelling, all 100 possible tracks generated for each de ployment using the 'probGLS' algorithm were in cluded as observed locations. Including all possible locations gives more weight to portions of the tracks where the location is more certain (because all possible locations are more clustered) and less weight to portions of the tracks where the location is less certain (because all possible locations are more dispersed). The mean standard deviation across iterations for any location estimate was 2.1° longitude and 2.0° latitude. Pseudo-absences were randomly sampled from ocean areas within 1000 km of any location collected within each tracking year. Pseudo-absences were sampled at a 1:1 ratio with observed locations (Barbet-Massin et al. 2012). Areas within 200 km of used locations collected within the same month were excluded from the selection area to ensure that pseudo-absences were outside of areas known to be occupied by murres at that time. For each tracking year, entire tracks from 70% of individuals were randomly assigned as training data and entire tracks from the remaining 30% of individuals were used as test data. Selection areas for pseudo-absences were determined separately for test and training data. Fig. S1 provides example maps of used locations and pseudo-absences.
We used random forests for our SDM, using the 'ranger' package in R (Wright & Ziegler 2017). The model was fit using the train function in the 'caret' package (Kuhn et al. 2021). We first ran hyperparameter tuning on a subsample of 5% of the data, considering combinations of the split rule ('gini' and 'extratrees'), minimum node size (5, 10, and 15), and 'mtry' (1−7). The area under the receiver operating characteristic curve (AUC) was used to identify the best combination of hyper-parameters on the sample data, and the selected hyper-parameters (mtry = 1, splitrule = extratrees, node size = 5) were used on the full model. The model was fit to the training data using repeated 4-fold cross-validation with 10 repeats, where each fold used 3 tracking years for model training and 1 tracking year for model testing.
Model accuracy was assessed using the AUC, F 1 scores (F1), and the Continuous Boyce Index (CBI). Model accuracy statistics (AUC and F1) are presented from cross-validation used in model fitting, which represents the accuracy in predicting probability of occurrence to unobserved years. We also present model accuracy for withheld test data, which represents accuracy in predicting to tracks of new individuals. AUC and F1 scores were calculated using the 'pROC' package (Robin et al. 2011), and CBI was calculated using the 'ecospat' package (Broennimann et al. 2021).
Variable importance measures were used to assess the relative contribution of each predictor to the model. Variable importance was calculated separately within each stage (moult, fall, winter, and spring), using the 'vip' package (Greenwell et al. 2020), to examine how habitat preferences changed through the nonbreeding period. Variable importance measures for each stage were scaled to values between 0 and 100. We calculated accumulated local effects (ALEs) to examine how each environmental variable influenced the predicted probability of use. The ALE shows the relative effect of each predictor variable on the model predictions, and this measure is not biased by correlation among predictor variables (Molnar 2019); positive values indicate an increase in mean probability of use and negative values indicate a decrease in the mean probability of use. ALEs were calculated using the 'ilm' package in R ( Molnar 2018); ALE values for each predictor variable are reported for each nonbreeding stage, in order to examine how habitat use changes among stages of the annual cycle.
In 2017 and 2018, tracked murres were included in a separate study examining the effects of increased reproductive investment on non-breeding behaviour. Prior to developing the SDM described above, we tested if treatments applied in that study had any influence on the SDM (see Text S1, Table S2). After confirming that there was no effect of treatment, all tracks from 2017 and 2018 were included in the final model.

Distributions by non-breeding stage
We compared predicted distributions for the 4 lifehistory stages: moult, fall migration, winter, and spring. Stage-specific distributions were calculated by predicting murre occurrence from the SDM at 3 d intervals over the period 1982−2019, then calculating the median predicted probability of use for each raster cell in each stage for each year. To quantify changes in the distribution of habitat, the stage-specific range areas were defined using the probability cut-off that included 90% of used locations. To quantify changes in the predicted distributions over time, we calculated a baseline range based on the mean distributions for the period 1982−1989. For each stage, we estimated 8 distribution measures: the total area, percentage overlap with the baseline distribution, median longitude, median latitude, western edge, eastern edge, northern edge, and southern edge. All distribution measurements were made using an Albers equal area projection with central meridian at 60°W and standard parallels at 45°N and 65°W. Range edges were calculated as the 5 th and 95 th percentiles of eastings (western and eastern edges) and northings (southern and northern edges) of all raster cells within the range. We used linear regression, with year as a predictor, to test for changes in distribution measures over time. Residual and q-q plots were used to check assumptions of the linear regression; Spearman's correlation tests were used to confirm linear regression results if normality assumptions were not met. A Bonferroni correction was used to account for multiple comparisons on the same seasonal distributions.
To investigate how changes in climate variables contributed to changes in stage-specific distributions, we compared mean values of ice cover, SST, air temperature, and wind speed between the 1980s (1982−1989) and the 2010s (2010−2019). We calculated the predicted distribution within each period, and identified regions where the predicted distribution declined, remained stable, or increased. We randomly sampled 50 points within these regions for each season, extracted the mean climate values for the 2 periods, and calculated the change in mean climate values from the 1980s to the 2010s. We used generalized least squares regression to determine how the change in climate varied among regions where the predicted range had declined, remained stable, or increased. Fixed variance weights for each region were used to account for unequal variance among regions. Residual and q-q plots were used to check assumptions of the regressions.

Fall and spring habitat phenology
Murres that breed in Hudson Bay migrate annually to the northwest Atlantic Ocean. To examine changes in the timing of habitat availability in Hudson Bay in spring and fall, we predicted the amount of suitable habitat in Hudson Bay at 3 d intervals for the years 1982 to 2019. Habitat area was quantified as the area of suitable habitat within the Hudson Bay marine ecoregion (Spalding et al. 2007); suitable habitat was defined as above, using the probability cut-off from model predictions that included 90% of used locations. We used a non-linear logistic regression curve to model the seasonal decline in habitat availability in fall and the increase in habitat available in spring as a function of DOY. Analysis was performed using the 'nls' function in R. We compared a null model with no effect of year to a model that included a trend with year.
All analysis was done using R version 4.1.0 (R Core Team 2021); p-values < 0.05 were considered significant for all parametric tests.

Tracking
We recovered data from 90 individuals during the non-breeding seasons over the 4 years of tracking (Table 1). Tracked murres followed a similar migration route and used the same wintering area during each year of tracking ( Fig. 1; Figs. S2−S5). Immediately following breeding, murres remained in Hudson Bay through the moult in September and October. Murres migrated through Hudson Strait to the Northern Labrador Shelf in November and December. During winter, murres spread out within the Labrador Sea, also reaching the Gulf of St. Lawrence, the Eastern Scotian Shelf, the East Greenland Shelf, and the Irminger Sea. They began migrating back through Hudson Strait and into Hudson Bay in April. The track for 1 murre was excluded during moult, because this individual stopped breeding in early August and migrated to the wintering area before moulting.
Mean fall and spring migration occurred on DOY 333 and 117, respectively. There was no difference in migration timing among the 4 years for fall (χ 2 = 2.15, p = 0.542) or spring (χ 2 = 1.39, p = 0.707). There was significant variation in the timing of migration among individuals each year (Fig. S6). Murres migrated over a period of 50 d in fall and 42 d in spring. Across the year, 95% of all tracked murres migrated between DOY 308 and 362 in fall and between DOY 90 and 141 in spring. We used these dates to summarize habitat use in 4 stages of the non-breeding period: moult (DOY 245−307), fall migration (DOY 308−362), winter (DOY 363−88), and spring migration (DOY 89−152).

SDMs
The SDM had good predictive performance across years and individuals. Cross-validation AUC was 98.7% and the F1 score was 94.2%. For the withheld test tracks from 30% of individuals, the AUC was 97.2%, the F1 score was 91.2%, and the CBI was 0.82. Ninety percent of used locations were located in areas with a predicted probability of use of 0.7 or higher; therefore, a probability of use of 0.7 was used as a cut-off for estimating suitable habitat.
SST and distance from colony were the 2 most important variables across all 4 stages (Fig. 2). In fall, winter, and spring, the probability of use was higher for SST < 4°C (Fig. 3); during moult there was very little effect of SST on probability of use. Probability of use declined with distance from the colony in all stages, with the strongest effects occurring during moult and spring. Bathymetry was an important predictor in all stages; shallow water (< 500 m) was preferred during moult, when murres are using shallow areas in Hudson Bay, and deeper water (> 2000 m) was preferred in fall, winter, and spring. Air temperature was an important predictor in fall, winter, and spring, with probability of use higher for temperatures < 4°C. In winter, probability of use declined with air temperatures less than −12°C. Probability of use de clined with increasing ice cover. Slope, DOY, and wind speed contributed the least to the model predictions.

Predicted distributions within each non-breeding stage
The non-breeding distribution of murres closely follows the receding ice, which 'pushes' murres out of Hudson Bay and Hudson Strait into the North Atlantic each year, and which murres then follow back into Hudson Bay each spring following melt ( Fig. 4 suitable non-breeding habitat from 1982 to 2019, we found significant changes in predicted distributions in fall and winter, with the greatest changes occurring in fall ( Fig. 5; Table S3). Overlap with the baseline range in fall declined by 0.7 ± 0.1% yr −1 (mean ± SE) ( Table S3). The fall range shifted west by 21.1 ± 4.1 km yr −1 , occupying all of Hudson Bay by the 2010s, and north by 5.0 ± 1.0 km yr −1 , with increased use of Hudson Strait, Foxe Basin, and Davis Strait (Fig. 4). Overlap in the winter distribution declined by 0.3 ± 0.08% yr −1 , and the winter distribution shifted west by 2.8 ± 0.7 km yr −1 (Table S3). The spring range shifted north by 3.8 ± 1.3 km yr −1 ; however, this trend was marginally non-significant (p = 0.06; Table S3). The most notable changes in spring distribution came from increased use of Hudson Strait and northern Hudson Bay. Results for all distribution measures are provided in Table S3. Ice cover declined, while SST and air temperature increased between the 1980s and the 2010s (Fig. 6; Table S4). The largest changes in ice cover occurred within the fall distribution, where the region of increasing suitable habitat had mean declines of 21.5 ± 0.51%. For fall, winter, and spring, ice cover declined more in regions of stable and increasing suitable habitat than in regions of declining suitable habitat. Regions of declining fall suitable habitat had the largest increases in SST (fall: 1.11 ± 0.03°C), while regions of increasing winter and spring suitable habitat had the smallest change in SST (winter: 0.27 ± 0.02°C; spring: 0.05 ± 0.01°C). The largest increases in air temperature occurred in areas of increasing suitable habitat during fall (4.83 ± 0.15°C). For fall, winter, and spring, regions with declining suitable habitat had less change in air temperature than regions with stable or increasing suitable habitat. For moult, air temperature increased more in regions of increasing suitable habitat than in areas with stable or declining suitable habitat. Overall changes in air temperature during moult and spring were of a smaller magnitude than during fall and winter. For winter and spring stages, wind speed tended to increase in areas of increasing suitable habitat and decrease in areas of declining suitable habitat.

Fall and spring habitat phenology
Fall habitat is available later than in the 1980s (F 910, 3 = 126.8, p < 0.001, Fig. 7). The date when suitable fall habitat reaches the midpoint of decline has increased by 0.38 ± 0.03 d yr −1 (t = 12.83, p < 0.001). There was no evidence that the amount of suitable habitat available at the start of fall (237 ± 1199 km 2 , t = 0.20, p = 0.84) or the rate of habitat decrease through fall (0.01 ± 0.02, t = 0.42, p = 0.67) have changed. According to this model, there were 200 000 km 2 of suitable habitat available until DOY 344 in 1982, whereas in 2019, the same amount of habitat was available until DOY 358 (mean change: 0.38 d yr −1 ).

DISCUSSION
Since 1982, the predicted distribution of murres from Coats Island during fall and winter has shifted north and west. SST, air temperature, and ice cover were important climate variables within our SDM, which accurately predicted the distribution of murres during the non-breeding period. Range expansion was associated with declining sea ice cover and warmer air temperatures, while range contraction was associated with increasing SST. The greatest changes in distribution have occurred in the fall, where habitat available in Hudson Bay has increased substantially. Other recent studies have predicted that unchecked anth ropogenic climate change will result in a northward shift in the winter distribution of multiple seabird species in the North Atlantic (Clairbaux et al. 2021); our study shows that climate change has already contributed to shifts in the non-breeding distribution of murres from Coats Island.
Our SDM approach assumes that the murre niche is at equilibrium and that we characterized the relevant components of the niche, and, thus, that the statistical relationship between environmental variables and murre distribution measured by our SDM can be extrapolated backward in time (Elith & Leathwick 2009). This latter assumption may not be valid if murres have more phenotypic flexibility than is captured in our training data, or if murres have adapted their habitat preferences in response to changing climate conditions. The SDM was developed using a limited set of climate predictors that are available for the entire period 1982−2019; these variables likely do not capture all elements of the niche of murres during the non-breeding period. In particular, nonbreeding distributions are likely driven by biotic interactions, like the distribution of prey, for which we did not have the relevant information (murre non- breeding season diet and the relevant fish and invertebrate distributions at depth are poorly known). Nonetheless, many fish and invertebrate distributions are strongly associated with SST and sea ice cover (Perry et al. 2005, Søreide et al. 2010; as such, a central assumption is that our model should capture these biotic interactions indirectly.
SDMs based on climate data are best suited for coarse-scale modelling of widely distributed, mobile species (Robinson et al. 2011). Finer-scale modelling of habitat use and ecological interactions by murres will require more precise tracking methods (e.g. GPS loggers), which are not yet feasible for year-round deployment on this species. However, despite the sea surface temperature uncertainty inherent in the GLS locations, we believe this SDM provides useful information on current and historic distributions of murres be cause the predicted distributions and habitat associations identified using this model agree with current knowledge of the species' biology (Moody & Hobson 2007, Gaston & Hipfner 2020. There have been substantial changes in the timing of modelled habitat availability for murres in Hudson Bay over the last 38 yr. The average date when fall habitat declines to less than 200 000 km 2 has in creased by 3.8 d decade −1 , and the date when spring habitat reaches 200 000 km 2 has advanced by 1.6 d decade −1 . Therefore, murres could spend 21 more days in Hudson Bay in 2019 than in 1982. The murres tracked in our study demonstrated significant among-individual variation in the timing of migration in and out of Hudson Bay in the fall (range: 50 d) and spring (range: 42 d). This indicates that individual condition or preference play an important role in the timing of migration for this species. Increased availability of habitat within Hudson Bay could have a positive effect on murres, by allowing greater flexibility in the timing of mi gration. Any potential benefit of in creased access to Hudson Bay, however, could be offset by changes in the food web associated with on-going changes to the marine climate (Hoover et al. 2013a,b) or increased competition (Piatt et al. 2020). Transient benefits of climate change and sea-ice loss have been documented for other ice-associated Arctic species (Laidre et al. 2020).
Studies of climate change impacts on phenology tend to focus on changes to spring phenology (Askeyev et al. 2007, Wolkovich et al. 2012, Parmesan et al. 2013, Gallinat et al. 2015, Kolářová & Adamík 2015. Our study highlights an example of an Arctic population that is ex periencing the greatest climate change induced impacts on phenology during the fall. Conditions during fall can play an important role in population demographics by influencing juvenile survival and determining condition of both adults and juveniles at the onset of winter, when many species face harsh environmental conditions. For murres, delayed fall migration could provide post-breeding adults with more time to complete their flightless moult and gain body reserves, while also providing juvenile murres more time to grow and gain experience flying and foraging before undertaking their first migration. It takes breeding murres approximately 50 d after egg-laying to raise a chick to nest departure (Gaston & Hipfner 2020) and an additional 35 d at sea before chicks are independent (Elliott et al. 2017). This limits the time during  Table S3. All distribution measurements were made using an Albers equal area projection with central meridian at 60°W and standard parallels at 45°N and 65°W which murres can begin nesting and successfully raise young. In creased time with suitable habitat within the breeding range could increase the window when murres can successfully breed. Our study supports the growing consensus that autumn phenology can be as sensitive as spring phenology to changing climates, especially for species which must undergo a feather moult before migration (Jenni & Kéry 2003, Brisson-Curadeau et al. 2020.
Many migratory species depend on matching the timing of breeding with seasonal peaks in resource availability in order to achieve successful breeding (Perrins 1970, Cushing 1990). The timing of ice-off is an important determinant of peak production in Arctic regions (Legendre et al. 1981); therefore, the changing spring conditions could influence the ideal timing of breeding. The modelled habitat conditions that correspond to migrating back to Hudson Bay are  occurring earlier now than 38 yr ago. Murres at Coats Island breed earlier in years with less ice cover in Hudson Bay (Gaston et al. 2005), and the mean laying date has advanced by 0.25 d yr −1 since 1990 (S. Whelan unpubl. data), indicating that murres are able to advance breeding in response to changing climate conditions. As murres now have 21 more days of suitable habitat available within Hudson Bay, and seem to be tracking the increased availability of spring habitat to avoid a mismatch, it may appear that climate change is advantageous. However, changing marine climate in Hudson Bay has altered prey composition during the chick-rearing period ) and may bring in competitors (Gaston & Woo 2008), leading to lower chick growth rates and ultimately fitness. Our model identified important non-breeding regions for murres from Coats Island, which show consistently high use at a multi-decadal scale, and could be important for marine spatial planning to mitigate the impacts of increased marine activities. During the non-breeding period, murres from Coats Island were most concentrated during moult when their range was restricted to central Hudson Bay. This is a critical stage in their annual cycle, when they are flightless and the males are caring for their dependent offspring. Moult is the time of year when distribution is determined more by the static variable distance from colony than by climatic variables, and also the period (outside of the breeding season) when, because of temporary flightlessness, murres have the least flexibility to alter their distribution in response to chang-ing conditions. This could make murres sensitive to changes in prey conditions at this time of year. Hudson Bay has relatively low levels of human activity, specifically shipping, which could pose risks to moulting murres through by-catch or oil pollution. However, shipping in Hudson Bay and Hudson Strait is expected to increase, including in regions that overlap with the moult distribution (Pizzolato et al. 2016, Dawson et al. 2018. The wintering habitat that we have identified for murres from Hudson Bay includes habitat for a significant portion of the global population of thickbilled murres (Frederiksen et al. 2016) as well as many other seabird species (Mallory et al. 2008, Frederiksen et al. 2012, Hedd et al. 2012, Linnebjerg et al. 2013, Fifield et al. 2017, Amélineau et al. 2018. Conditions during winter probably play an important role in population regulation of murres (Gaston 2003). Across their annual range, the wintering area is the region where murres from Coats Island experience the greatest overlap with anthropogenic threats such as shipping, fisheries by-catch (Davoren 2007), chronic and acute oil pollution (Wiese & Ryan 2003), and hunting (Gaston & Robertson 2010). Modelling winter habitat use, as we have done here, is an important step in developing marine spatial planning of offshore wintering areas for seabirds. The extent to which murres can shift their winter distribution northward may be affected by day length. As murres forage primarily, though not exclusively, during daylight (Regular et al. 2011, Dunn et al. 2020, northern range shifts would involve a decrease in the amount Day of year Spring habitat (1000 km 2 ) of time available for foraging in daylight. A comparison of winter energy expenditure and dive behaviour of thick-billed murres and common murres wintering in regions with and without polar night showed that wintering in an area with polar night resulted in higher daily energy expenditure and reduced foraging opportunities for common murres .
Diversity of responses among populations within a species can enhance its resilience (Sydeman et al. 2015). Our analysis focussed on habitat changes for a population of murres close to the southern edge of their range. Other populations in northern parts of the range faced with similar environmental changes may experience a different response. For example, colonies in the Canadian High Arctic are more constrained by sea ice during the breeding season (Gaston et al. 2005). Earlier ice-off dates around these colonies could lead to improved reproductive success for these populations. Negative effects for populations within some portions of the breeding range may be offset at the species level by positive effects in other regions. However, to the extent that large portions of the murre population share common wintering areas, and conditions during the winter play an important role in adult survival (Frederiksen et al. 2016), this could create greater sensitivity to climate change at the species level.
Mapping the distribution of populations outside of the breeding season is a key priority for the conservation of marine birds. Long-term data on winter distributions of pelagic seabirds were, until recently, limited to data collected from band recoveries and vessel surveys, which are biased in the Arctic by low spatial and temporal coverage. In contrast, our SDM was effective at predicting stage-specific distributions of murres across years and individuals and is an improvement over the more common approach using utilization distributions to map species' ranges, which are limited by the relatively small number of tracked individuals. The SDM approach is more likely to identify the entire distribution, which could otherwise be missed. In dynamic marine environments, an SDM can account for significant inter-and intra-annual variation in habitat by predicting distributions over multiple time periods. Our study demonstrates the advantage of including year-round tracking as part of long-term monitoring programmes to facilitate improved understanding of non-breeding distributions, habitat requirements, and the effects of environmental variability and climate change on population demographics (Carneiro et al. 2020).