Unexpected sources of uncertainty in projecting habitat shifts for Arctic shorebirds under climate change

The rapidly changing Arctic is ideal for investigating uncertainties in climate projections. Despite the challenges of collecting data in this region, an unprecedented large‐scale survey of shorebirds has been conducted over the last 30 years. Our study aimed to (1) develop probabilistic estimates for the change in suitable habitat for 10 Arctic shorebird species in Canada by 2075 and (2) assess the contribution of modelling decisions to the uncertainty in these estimates.


| INTRODUC TI ON
The location of many species' habitat has been shifting in response to climate change, and forecasting these shifts has been an active area of research for more than 25 years (Carey, 1996;Peterson et al., 2002).One of the most common ways to do this is with correlative approaches.Habitat is first defined by relating contemporary species occurrences to the resources and conditions in their environment, and then these relationships are used to project future habitats.Projections of how species' habitat will shift in response to climate change have been useful for identifying the scale of change in species ranges that may occur in response to climate change (e.g.Loarie et al., 2008;Morato et al., 2020).They are also useful for identifying which species, habitats, and regions may be most sensitive to climate change (e.g.Chambault et al., 2022;Dyderski et al., 2018).
However, these projections are hampered by a high degree of uncertainty from numerous sources (Peterson et al., 2018).All habitat suitability models are affected by uncertainty from sources such as limited sample sizes, measurement error, incomplete detection, imperfect understanding of what makes habitat suitable, lack of data for important predictor variables, and uncertainty in which model structures best represent the shape of ecological relationships (Barry & Elith, 2006;Beale & Lennon, 2012).However, for future projections of habitat suitability, additional sources of uncertainty must also be considered.These include uncertainty in how human behaviour will shape carbon emissions and land use, how the earth's climate system is modelled, how climate changes will affect habitat, how species-habitat relationships may shift over time, and an inability to directly evaluate model accuracy using independent data (Heikkinen et al., 2006).There are many ways to categorize these sources of uncertainty, but two broadly important categories are epistemic uncertainty, i.e. uncertainty due to lack of knowledge, and stochastic uncertainty, i.e. uncertainty due to natural variability (Kirchner et al., 2021;Regan et al., 2002).
A substantial body of scientific literature has developed to address epistemic uncertainty by identifying how modelling decisions affect future projections of habitat suitability.These papers have explored the relative importance of sources of uncertainty from model structure, such as modelling algorithms (Diniz-Filho et al., 2009;Pearson et al., 2006), approaches for variable selection, model parameter settings (Peterson et al., 2018), thresholding procedures (Steen et al., 2017), and choice of predictor variables, which have been limited to various permutations of climate variables for philosophical or practical reasons (Braunisch et al., 2013;Synes & Osborne, 2011;van de Pol et al., 2016, but see Stralberg et al., 2015).
Other sources of model uncertainty arise from decisions about the data inputs used to create projections, such as the global circulation models that describe spatiotemporal patterns in climate (GCMs; Braunisch et al., 2013) and the carbon emissions scenarios that force these models (Beaumont et al., 2008).Biotic interactions have been argued to be important for predicting how species distributions will respond to climate change (Van der Putten et al., 2010).They have historically been left out of future projections of species habitat, but they are gaining popularity as suitable data becomes available (e.g.Norberg et al., 2019;Ovaskainen et al., 2017).
In this paper, we compare sources of uncertainty in projections of how climate change may affect the breeding habitat of Arcticbreeding shorebirds.This system is ideally suited for exploring the sources of uncertainty in climate-related projections.Uncertainty is particularly important in the Arctic, as many sources of uncertainty are exaggerated here.Expected climatic changes are more extreme than in other regions; temperatures in the Arctic are rising three times faster than the global average, and impacts are expected to be rapid and transformational (AMAP, 2021).At the same time, there have been fewer and less robust studies of the effects of climate change here relative to other terrestrial biomes, as consistent data collection efforts are much more limited, both for wildlife observations and for habitat predictors such as climate and land cover variables (Aronsson et al., 2021).However, in the last 30 years, shorebirds have been well monitored in an unprecedented Arcticwide survey (Bart & Johnston, 2012).
These Arctic-breeding shorebirds appear to be shifting their ranges northwards in response to recent climate change (Anderson, Fahrig, Rausch, Martin, et al., 2023).A qualitative risk assessment of the vulnerability of Arctic-breeding shorebirds to climate change in North America predicted that many of these species will lose more than half of their suitable breeding habitat, at a high level of subjective confidence (>70%, Hector Galbraith et al., 2014).A correlative model projecting circumpolar changes in Arctic shorebird breeding habitat suggested that, globally, 66-83% of shorebird species could lose half of their suitable breeding habitat by 2070 (Wauchope et al., 2017).However, neither of these studies provides a comprehensive, quantitative picture of the uncertainty around these future projections of Arctic shorebird breeding habitat, which we hypothesized should be very high.
Our study had two objectives.The first was to develop probabilistic estimates of how the area of suitable habitat for 10 species of Arctic-breeding shorebirds in Canada will change by 2075.The stochastic uncertainty from natural variation in ecological systems is difficult to predict; therefore, one of the best ways to represent this uncertainty is as a probability distribution (Elith et al., 2002).The second objective was to assess how modelling decisions contribute to the uncertainty in these estimates, assessing a range of model configurations that are currently used in habitat suitability studies.
Our assessment of uncertainty included six distinct classes of modelling decisions, three of which, to our knowledge, are not well studied in papers about uncertainty in species' responses to climate change.
First, we determined how the method for selecting habitat variables in habitat suitability models influences uncertainty, comparing a full model and two variable selection methods based on cross-validation (see Dormann et al., 2008 for uncertainty from AIC-based methods).

While older stepwise regression techniques have been out of favour
for some time (Whittingham et al. 2006), methods based on spatial cross-validation are currently being debated as a modern alternative (Meyer et al., 2018;Wadoux et al., 2021).There is also debate on whether it is appropriate to avoid variable selection, instead using a "full-model" approach, including a large number of highly correlated variables (Beale & Lennon, 2012;Smith et al., 2009).Second, we compared the consequences of the pool of candidate variables available for selection.We compare such climate-only models and models including additional non-climate variables, which are important predictors of Arctic shorebird distributions (Anderson, Fahrig, Rausch, & Smith, 2023).Third, we tested how assumed dispersal limits in projections of future vegetation affect uncertainty.In addition to these three classes of modelling decisions which have not been well studied, we explored three more commonly tested classes of modelling decisions: the modelling algorithm, carbon emissions scenario, and global circulation model.

| Overview
Our study had two goals: (1) to develop probabilistic estimates of how the area of suitable habitat for 10 species of Arctic-breeding shorebirds in Canada will change by 2075, and (2) to assess how modelling decisions contribute to the uncertainty in these estimates.We varied 6 different modelling decisions, giving 216 unique projections of future habitat for each species.For each projection, we calculated the proportional change in suitable habitat as the future suitable area divided by the present suitable area.We used a bootstrapping approach, creating a probability distribution for the proportional change in suitable habitat estimated for each unique projection.We then combined these to create an overall probability distribution for each species that captures both within-and between-model uncertainty.We compared sums of squares to determine the relative contributions of different modelling decisions to uncertainty in our estimates of the proportional change in suitable habitat.The analyses are summarized in Figure 1 and described below, as well as in the ODMAP (Overview, Data, Model, Assessment and Prediction;Zurell et al., 2020) document in Appendix S3.All analyses were performed in R version 4.2.3, using packages tidyverse 2.0.0, caret 6.0-94, and terra 1.7-29.See Anderson et al. (2024) for the full dataset used in our analysis.

| Shorebird surveys
Plots were surveyed for shorebirds across the Canadian Arctic as part of the Arctic Program for Regional and International Shorebird Monitoring (PRISM; Bart & Johnston, 2012), which covers all of Arctic North America as defined by the Circumpolar Arctic Vegetation Map (Walker, 2000).The Canadian PRISM data used here includes 1854 plots across the Canadian Arctic (Figure S1.1), an area of 3.5 million km 2 .Given this massive area, surveys were completed over a number of years, from 1994 to 2018, with a majority of the sampling effort concentrated between 2003 and 2018.The Canadian Arctic was divided into 12 regions based on physiognomy and logistical considerations.Within each region, plots were selected by random sampling, stratified by habitat.The proportion of plots in wet, moist and dry habitats varied by region depending on habitat abundance (See Bart & Johnston, 2012 for more details), but on average 16%, 29%, and 56% of plots were in wet, moist, and dry habitats respectively.In all years, surveys were conducted between June 18 and July 15.This corresponds to the latter part of shorebird courtship and the early part of incubation, when detectability was highest, as breeding territories could be identified from the birds' territorial displays.Although we use the uncorrected densities in our analyses, the PRISM survey also includes a subset of plots that are surveyed intensively, so that shorebird detectability could be estimated using double-sampling methods.The estimated detectability rate is 1.13, with no significant difference between species, indicating that occasionally birds visiting the plot were misidentified as breeding (J.Bart and P. Smith, unpublished).Plots were typically 12-16 ha (400 m × 300-400 m).We excluded 23 surveys that were conducted in 1996 when extensive flooding led to very atypical breeding conditions.We also excluded 36 plots for which observers recorded nearby human development, which is rare across the study area and could potentially alter shorebird habitat preferences.Finally, 480 plots were excluded because not all necessary environmental predictors were available for them (for example, plots in areas close to shore which were in locations classified as ocean in the WorldClim dataset).Our final dataset consisted of surveys from 1854 plots.For 145 plots that were visited in more than 1 year, we randomly selected 1 year of data to include in the final dataset.
Plots were surveyed for breeding birds following PRISM protocols (Bart & Johnston, 2012).Two surveyors searched the plot walking straight-line transects, covering a breadth of 50 m with the observers situated 25 m apart, using a GPS to ensure complete coverage of the entire plot.The survey of a single plot took approximately 90 min.Observers recorded the number and species of all birds observed within each plot.Altogether we observed 19 species of shorebirds, 10 of which had data sufficient to minimize overfitting in the models of suitable habitat (observed in at least 30 plots and containing at least 60% of their North American range within the study area): American Golden-Plover (Pluvialus dominica), Baird's

| Environmental predictors
For our habitat suitability models, we used two categories of environmental predictors of shorebird occupancy: climate predictors found in 'climate-only models', and additional predictors in Additional predictors in the 'climate+additional' models of shorebird occupancy included topography, geology, and land cover variables.
The topography predictors we included were elevation, the standard deviation of elevation, and distance to the coast.Elevation data were from WorldClim 2 (1 km resolution), from which we also derived the standard deviation of elevation.We calculated the distance to coast from coastline data from Natural Earth (Patterson & Vaughn Kelso, 2022).To represent geology, we used substrate chemistry data from the Circumpolar Arctic Vegetation map (CAVM; Walker et al., 2005), which includes 3 categories: acidic (pH < 5.5), circumneutral (pH 5.5-7.2) and carbonate (pH > 7.2).We used land cover data from the Circumpolar Arctic Vegetation Map (CAVM; Walker et al., 2005), which derived 17 land cover classes from AVHRR data from 1993 to 1995.We used 13 of these land cover classes, including barren, graminoid tundra, prostrate-shrub tundra, erect-shrub tundra, and wetland landcover classes.We excluded areas classified as glacier, mountain, water, and non-Arctic areas from our analyses, as these areas were not surveyed and are not suitable for shorebird habitats.Both CAVM datasets are vector maps with a minimum polygon edge size of 14 km.
For future projections of shorebird habitat, we used the above predictors for topography and geology, as well as future projections of bioclimatic variables and land cover for 2040-2060.We used several different sets of bioclimatic variables from the WorldClim 2 dataset (more details below).Future projections of the CAVM land cover maps described above were created by Pearson et al. (2013), based on climate, geological substrate, and limitations on the distance that the tree line can shift northwards.Pearson et al. (2013) did not make future projections of land cover classes that do not have enough soil to support vegetation (barren land cover classes), because it is expected that these areas will remain unchanged.They also did not make projections for land cover classes that are highly influenced by hydrological processes (wetland land cover classes), because the changes could not be adequately predicted using data available at the time.For these areas, we based both recent and future habitats on the original CAVM classification.All recent and future predictors were resampled to a 5 km raster grid using bilinear interpolation for continuous data and nearest neighbour resampling for categorical data.
Future projections created using habitat suitability models may perform poorly when extrapolating to environmental conditions which were not included in the training data (Fitzpatrick & Hargrove, 2009).To minimize extrapolation in our future projections, we mapped the range of values for both the recent and future environmental predictors to look for conditions that were far outside the sampled range of values.Future precipitation regimes in the Arctic Cordillera, and the eastern portion of Baffin Island, are predicted to diverge considerably from the climate conditions observed in areas that were sampled in this study.We therefore excluded these areas from predictions.

| Modelling decisions
To account for and assess how modelling decisions contributed to uncertainty in our estimates of the proportional change in suitable shorebird habitat, we varied six different modelling decisions in our analyses.Three of these decisions concerned the structure of our habitat suitability models, with 12 possible configurations, and three of these decisions concerned the environmental input scenarios for future projections, with 18 possible configurations.This created a total of 216 unique modelling decision configurations for each species (Figure 1).
The three modelling decisions about model structure that we varied were the pool of candidate variables, the modelling algorithm, and the selection of variables.For these three decisions, we evaluated options intended to represent a wide spectrum of decisions that are currently made in habitat suitability studies.For the pool of potential variables, we used with either 'climate-only' candidate predictor variables (e.g.Dyderski et al., 2018) or 'climate+additional' candidate predictor variables (e.g.Heikkinen et al., 2007; see Section 2.3).For the modelling algorithm, we created habitat suitability models with either generalized linear models (GLM) with a binomial distribution (e.g.Chiaverini et al., 2023), or a Random Forest F I G U R E 1 Methodology used to estimate the probability distribution of the proportional change from 2010 to 2075 in suitable breeding habitat for each of 10 shorebird species.In Step 1, we built 12 habitat suitability model structures for each species, with all possible configurations of 3 modelling decisions (the pool of candidate variables, the modelling algorithm, and the selection of variables).
In Step 2, we compiled 18 environmental input scenarios, with all possible configurations of 3 modelling decisions (GCM, carbon emissions scenario, and tree line dispersal).In Step 3, we created present and future projections for each of the 216 model structure × environmental input configurations.We then calculated the proportional change in suitable habitat as the future suitable area divided by the present suitable area.In Step 4, we repeated Step 3 with 500 bootstrapped resamples of the shorebird occupancy data for each of the 216 model configurations, estimating a probability distribution for the proportional change in suitable habitat.We combined the 216 probability distributions to create a single probability distribution for each species.
algorithm (Cushman & Wasserman, 2018).For selection of variables, beginning with the pool of candidate predictors, we used three different methods for selecting which variables to include in the final model: a full model including all candidate variables (e.g.Hof et al., 2012;Langham et al., 2015;Pearson et al., 2013;Wauchope et al., 2017), a reduced model using forward variable selection with random cross-validation (e.g.Hao et al., 2019), or a reduced model using forward variable selection with spatial cross-validation (e.g.Chambault et al., 2022).For both reduced models, we used the R package CAST forward feature selection method to select the most important variables by testing all possible two-pair permutations of variables (Meyer et al., 2023).We used Cohen's Kappa statistic as our assessment metric for model fit and accuracy, which measures the extent to which the agreement between observed and predicted values is higher than that expected by chance alone (Liu et al., 2011).A Kappa statistic of 0 indicates a level of agreement expected by chance, while a Kappa statistic of 1 indicates perfect agreement between observed and predicted values.The model with the highest Kappa statistic was selected, and the process was iteratively repeated, adding each of the remaining variables, until none of the remaining variables increased model performance.The Kappa statistic was calculated using 4-fold cross-validation in one of two ways: in random cross-validation, the data were randomly split into 4 folds, while in spatial cross-validation, the data were split into 4 folds based on the 12 survey regions, previously described in the Shorebird Surveys section (more details in Meyer et al., 2018).We excluded species from our analysis that had low sample sizes (<30 presences) and Kappa close to 0.
The three modelling decisions about model inputs that we varied for future projections were the global circulation model (GCM), the carbon emissions scenario, and the maximum distance of tree line dispersal.For these three decisions, we evaluated the options as close as possible to those presented in Pearson et al. (2013) to ensure that we had land cover predictions to match each configuration of the modelling decisions.Thus, following Pearson et al. (2013), the GCMs we used were CanESM5, HadGEM3-GC31-LL, and ACCESS-CM2 (IPCC, 2022).Pearson et al. (2013)

| Projecting change in suitable habitat
We created estimates of the area of suitable habitat for 10 Arcticbreeding shorebird species for two time periods, the recent past and the future.First, we estimate the area of suitable habitat for the period when the shorebird observation data were collected, 1994 to 2018.We refer to this time period as the estimates for 2010, which is the median year for the observations.In our models, there are lags between the observation period for shorebirds, and the observation period for the climate data  and the land cover data (1993)(1994)(1995).We chose these data sets because they are the only options available that would allow us to include future land cover predictors and the influence of tree dispersal on shorebirds, which we identified as being important variables in Anderson, Fahrig, Rausch, and Smith (2023).We also believe there is likely a true lag in the responses between climate, landcover, and shorebird occupancy, so these time lags are a reasonable modelling assumption.The future climate predictions we use are for 2040-2060, therefore we carry forward the 25-year time lag, and estimate the area of suitable shorebird habitat for 2065-2085.We refer to this time period as the estimates for 2075.
Probabilistic methods are a comprehensive and effective way to deal with compounded effects of uncertainty from within models, among candidate models, and among different future climate scenarios (Wenger et al., 2013).Our process for projecting the proportional change in suitable habitat for each shorebird species involved 4 steps (Figure 1).In Step 1, we built 12 habitat suitability model structures for each species (see Modelling Decisions).In Step 2, we compiled 18 environmental input scenarios (see Modelling Decisions).In Step 3, we created present and future projections for each of the 216 model structure x environmental-input configurations.We then calculated the proportional change in suitable habitat as the future suitable area divided by the present suitable area.The probability threshold for classifying a grid cell as suitable or unsuitable was calculated as the threshold that maximized the Kappa statistic.In Step 4, we repeated Step 3 with 500 bootstrapped resamples of the shorebird occupancy data for each of the 216 model configurations, estimating a probability distribution for the proportional change in suitable habitat.We combined the 216 probability distributions to create a single probability distribution for each species.Measurements of proportional change follow a non-linear scale (i.e.half as much habitat is a proportional change from 1 to 0.5, while twice as much habitat is a proportional change from 1 to 2).Therefore, in our models, we compared the logarithm of predicted proportional change so that increases and decreases would have values of comparable magnitudes (i.e. for half as much habitat, log(0.5)= −0.69,and for twice as much habitat, log(2) = 0.69).

| Estimating uncertainty due to modelling decisions
To estimate how much of the uncertainty in the proportional change in suitable habitat was related to modelling decisions, we used their proportional contributions to the sums of squares in two ANOVA models.In both ANOVA models, the proportional change in suitable habitat was the response variable.In the first ANOVA, we included two fixed effects: the overall influence of modelling decisions as a factor with 216 unique levels reflecting all possible configurations of model structure and environmental input data, and species to control for variation among species.In the second ANOVA, we included each of the six classes of modelling decisions as a fixed effect, to compare their relative effects on prediction uncertainty, and we included species as a random effect.We also included an interaction effect between the modelling algorithm and the selection of variables, as modelling algorithms can differ in their sensitivity to collinear predictors in a full model approach (Dormann et al., 2013).
We also present a set of ANOVA models designed to look at how modelling decisions influenced the likelihood of large changes in suitable habitat, which may be a more appropriate way to describe the results of highly uncertain future projections.In these ANOVA models, the response variable was a binary variable predicting whether the species will have less than half as much suitable habitat in the future, or twice as much suitable habitat in the future.Again, we included each of the 6 modelling decisions as a fixed effect, and species as a random effect.

| Consensus maps of suitable habitat
For each of the 216 model configurations for each species, we created maps of suitable habitat for 2010 and 2075.For each time period, we created a consensus map, illustrating the proportion of the 216 models that predicted a grid cell to be suitable habitat, with suitable habitat identified as defined above, using the probability threshold that maximized the Kappa statistic.With the present and future consensus maps for each species, we assessed which areas were consistently predicted to be suitable in each time period, how these areas shifted over time, and which regions were consistently predicted to be suitable habitat for multiple species.

| RE SULTS
Compared to the area of suitable habitat predicted for 2010, the median predicted area of suitable habitat for 2075 ranged from 15% for Baird's Sandpiper to 178% for Semipalmated Sandpiper (Figure 2).
Across the wide range of modelling decisions considered here, 8 of 10 species were more likely to have less suitable habitat in 2075 than in 2010 (Baird's Sandpiper, Pectoral Sandpiper, White-rumped Sandpiper, American Golden-plover, Buff-breasted Sandpiper, Red Phalarope, Black-bellied Plover, and Stilt Sandpiper; Table 1).These same eight species had at least a 30% chance of losing half their suitable habitat by 2075.Two species were more likely to have more suitable habitat in 2075 than in 2010 (Dunlin, and Semipalmated Sandpiper; Table 1).These species, as well as Stilt Sandpiper, had at least a 30% chance of having twice as much suitable habitat in 2075.
Although the weight of evidence clearly suggests a loss of habitat is likely for the majority of these species, the uncertainty is substantial.The mean Kappa statistic across species was 0.20 (0.12-0.35), indicating meagre predictive power (Table 1, Table S1.1).
For four species, there was at least a 20% chance of having both half or twice as much suitable habitat in 2075 (Pectoral Sandpiper, Buff-breasted Sandpiper, Black-Bellied Plover and Stilt Sandpiper; Table 1).However, the degree of uncertainty varied considerably between species.Buff-breasted Sandpiper had the widest probability distribution (SD of 4.86) was five times broader than that for Semipalmated Sandpiper (SD of 0.88; Figure 2).The standard deviation in the proportion of habitat loss was higher for species that were rarer in our samples (β = −0.01,t value = −2.83,p = .02).
In our ANOVA model testing, the overall influence of modelling decisions on proportional change in suitable shorebird habitat, we found that most of the uncertainty in our estimates of the proportional change in suitable habitat was unexplained by the factors we tested, and only 10% of the total sum of squares was related to modelling decisions (Figure 3a).However, it is noteworthy that variation we observed among modelling decisions is comparable to the variation we observed among species (12% of the total sum of squares), so it is also not inconsequential.See Discussion for a review of some of the sources of uncertainty that were not included in this model.
In our ANOVA model of how the six individual modelling decisions affect predictions of proportional change, the largest source of uncertainty was related to the choice of approach for selecting the variables that were retained in the habitat suitability models, accounting for 67% of the sum of squares (Figure 3b).Choice of GCM

TA B L E 1
The probability that the area of suitable habitat for each shorebird species in 2075 will decrease compared to 2010 (proportional change <1), will be less than half of the area of 2010 (proportional change < ½), or will be more than twice the area in 2010 (proportional change >2).For each species, these values represent the cumulative probabilities beyond these thresholds, illustrated in Figure 2 as the areas above and below the red lines.accounted for 14% of the total sum of squares.The remaining modelling decisions contributed a small amount of uncertainty to the sum of squares: the choice of carbon emission scenario at 2%, and the choice of modelling algorithm, the pool of candidate variables at 1%, and the maximum distance of tree line dispersal at less than 1%.
The interaction between the choices of variable selection methods and modelling algorithms accounted for 15% of the uncertainty.
How we selected the variables that were retained in the habitat suitability models also had the by far strongest influence on the uncertainty in predicting whether a species will lose or gain a large amount of habitat, respectively, accounting for 79% and 77% of the explained sum of squares in these models (Figure 3c,d).The choice of GCM was the second largest source of uncertainty our model of the probability that species would lose half of their suitable habitat, at 10% (Figure 3c).The choice of modelling algorithm was the second largest source of uncertainty in our model of the probability that species would have twice the amount of suitable habitat, at 17% (Figure 3d).The other modelling decisions accounted for less than 3% of the sum of squares in both models.The interaction between the choices of variable selection methods and modelling algorithms was less important here; it accounted for less than 5% of the uncertainty in both the models of habitat loss and gain.
By plotting the probabilities of predicted half or twice as much suitable habitat associated with the choices for each modelling decisions across all 10 species, we can see that there was substantial overlap in the predicted probabilities for most modelling decisions (Figure 4).The largest difference in predicted probabilities is when different methods were used to select included variables.Both methods selected an average of four climate predictors, and in the 'climate+additional models', an average of one additional predictor (Table S2.1), compared to 19 climate variables and potentially 5 additional variables when all variables were included.Species are predicted to lose much less, or gain much more habitat when all candidate variables are included, compared to either of the two variable selection approaches was used (Figure 4).Likewise, we can see that GLMs consistently estimated more optimistic outcomes than Random Forest models.
Interestingly, the uncertainty arising from the carbon emissions scenarios was only 2% or less of the sum of squares in our models of how individual modelling decisions affected the proportional change in suitable habitat, and the probability of predicting half or twice as much suitable habitat (Figure 3b-d).Compared to the SSP 245 scenario ("middle-of-the-road"), the proportional change in the area of suitable habitat in the SSP 585 scenario ("fossil-fueled development") ranged from 4% more habitat to 13% less habitat, depending on the species (mean = 6% less habitat).This translated to an average 4% higher probability of predicting half the area of suitable habitat in the SSP 585 scenario (0-10% by species), and an average 1% higher probability of predicting twice the area of suitable habitat (−1 to 2% by species).For context, the mean average annual temperature in our study area was −12.6°C from 1970 to 2000.By 2040-2060, it is predicted to increase to −6.5°C in the SSP245 scenario, and −5.3°C in the SSP 585 scenario.
Mapping the consensus in the projections of suitable habitat across all model configurations illustrates a general trend where the location of suitable habitat for most species will shift north and west by 2075 (Figure 5).Even for species such as Stilt Sandpiper, Blackbellied Plover, and Red Phalarope, which were predicted to be less likely to lose or gain a large amount of habitat, the location of their suitable habitat was predicted to shift considerably (Figure S1.2).

| DISCUSS ION
Our results showed that many Arctic-breeding shorebird species have a considerable risk of losing over half of their suitable breeding habitat in Canada, but this projection is much less certain than has been described previously.The probability distributions we projected for these 10 Arctic shorebird species showed that decreasing area of suitable habitat is the most likely prospect for 8 of the species.However, for most species, a wide range of scenarios from major losses to modest changes to major gains in suitable habitat were all possible.Previous studies of shorebird habitat shifts do not show this high degree of uncertainty.Wauchope et al. (2017) predicted that globally, 66-83% of shorebird species would lose at least half their area of suitable habitat by 2070.They combined 3 GCMs into a single consensus model, and the only uncertainty addressed in this paper is a range of values creates from 2 emissions scenarios.
We found both GCM and emissions scenarios to be relatively minor This study also included Alaska, where there is no land further north to accommodate range shifts northwards.
There are a number of reasons why Arctic-breeding shorebirds are likely to be vulnerable to climate change: they depend on ecological synchronicities with short peaks of invertebrate prey availability; they migrate very long distances; and they specialize on particular habitat types (Galbraith et al., 2014).Tundra habitats are experiencing a faster rate of climate change than those at lower latitudes, due to positive feedbacks from changes in snow and ice dynamics (Serreze & Barry, 2011), leading to increased shrub size and cover (Myers-Smith et al., 2011), which is reducing the habitat suitability for most shorebird species (Saalfeld et al., 2013).
We were surprised that, of the modelling decisions we assessed, the largest source of uncertainty in projections of habitat change was from our choice of method for selecting the variables that were retained in the habitat suitability models.Specifically, when all candidate variables were included, species were predicted to lose much less, or gain much more, habitat than when a variable selection approach was used.This was an unexpected finding, as compared AIC-based methods.We found that variable selection was a much larger source of uncertainty than either the modelling algorithm or the GCM, modelling decisions that are regularly identified as significant sources of uncertainty in other studies (Buisson et al., 2010;Diniz-Filho et al., 2009;Naujokaitis-Lewis et al., 2013;Pearson et al., 2006).Our results showed that uncertainty stemming from variable selection depends on the modelling algorithm used.The full model approach may be more applicable for machinelearning methods, which may better cope with high-dimensional, collinear datasets (Kampichler & Sierdsema, 2018), but these methods may be less appropriate for GLMs.The method used for estimating relative effects can also influence the results; for example, variable partitioning methods such as we used here may be more strongly biased by correlations between predictors in 'full model' approaches than other methods for assessing relative effects (Smith et al., 2009).The additional uncertainty introduced by collinear variables in the full model approach is one area where uncertainty in future projections can be reduced by following best practices in structuring habitat suitability models (De Marco & Nóbrega, 2018;Sofaer et al., 2019).Forward selection approaches, which are much more conservative that a full model, purported to only include variables with actual predictive power beyond the training locations (Meyer et al., 2019).There is no consensus on the best way to select which variables are included in a model.Both full model and stepwise selection approaches have been proposed to avoid overfitting models to the samples data (Harrell, 2001;Meyer et al., 2018).We recommend that future uncertainty analysis take variable selection into account, and that further research should quantify uncertainty in habitat suitability projections across a larger array of variable selection techniques.
The limited influence of the choice of carbon emissions scenarios on the results is also notable.Similar to many other studies of uncertainty in projections of species distribution and habitat in response to climate change, we found that uncertainty from other modelling decisions greatly exceeds uncertainty from carbon emissions scenarios, the only modelling decision that reflects real-world decisions that we have some control over rather than a consequence of climate change or strictly a methodological uncertainty (Buisson et al., 2010;Garcia et al., 2012;Peterson et al., 2018;Synes & Osborne, 2011).The small differences in the probability of species predicted to lose or gain large quantities of suitable habitat with different carbon emissions levels suggests that many of the harmful effects of climate change on wildlife distributions are difficult to halt or reverse given plausible emissions scenarios.
It is important to note that we found that uncertainty in our projections stemming from modelling decisions was relatively small compared to uncertainty from other sources.Our study, like all correlational projections, is constrained by current observed species-habitat relationships.We found that there were higher levels of uncertainty for rarer uncommon shorebirds.One limitation of habitat suitability models is that for rare species, many areas of suitable habitat are likely to be unoccupied (Doncaster et al., 1996).Another limitation is that occurrence data may not adequately capture the niche of a species, either due small sample sizes or biases in sampling (Loiselle et al., 2007).
Opportunistic species, such as Pectoral Sandpipers, have low site fidelity, and selection of breeding sites may be more closely tied to annual conditions rather than long-term climate averages (Saalfeld & Lanctot, 2015).Habitat suitability models are also limited in their ability to capture some important ecological processes that operate at a range of spatial scales to influence whether suitable habitat is actually F I G U R E 5 Consensus maps of the areas predicted to be suitable habitat for 10 Arctic shorebird species in 2010 and 2075.The value in each cell is the number of species for which at least 50% of the 216 model configurations of 6 modelling decisions indicated suitable habitat.
used by a species, such as biotic interactions, dispersal limitations, and disturbance (Guisan & Thuiller, 2005).There is still much we do not know about the current distribution of shorebirds, and although we have some understanding of what defines suitable habitat for shorebirds, our ability to predict their distribution has limitations (Anderson, Fahrig, Rausch, & Smith, 2023;Cunningham et al., 2016;Saalfeld et al., 2013).Any limitations in models of present habitat suitability persist in the future projections they are used to create.Similar to Hegyi and Laczi (2015), we found that the choice of variable selection method and other statistical modelling decisions have a smaller impact on the results of ecological models than the quality of the data.
Using probabilistic methods was key for understanding the relative importance of modelling decisions as a source of uncertainty.
Araujo and New (2007) describe probabilistic forecasting as the 'end game' of multi-model ensemble forecasting, but these methods are still surprisingly infrequently used in projections of how species' habitat will shift in response to climate change.By explicitly describing uncertainty in future projections of species habitat explicit, it can be incorporated into conservation planning to make decisions that are robust to uncertainty (e.g.Bagchi et al., 2013;Carvalho et al., 2011;Fuller et al., 2008;Sofaer et al., 2019).For example, if there is consensus between multiple modelling methods that a particular site has a high conservation value, its protection can be prioritized over sites with greater uncertainty.
We have illustrated a high degree of uncertainty in how suitable breeding habitat for Arctic-breeding shorebirds will change under scenarios of climate change, but in fact the true uncertainty is even greater when we consider factors that were not possible to include in our models.For example, the future land cover predictions we included did not address changes in hydrology (Pearson et al., 2013).
Shorebirds are likely to be sensitive to such changes, as many species rely on moist, lowland habitats (Cunningham et al., 2016;Saalfeld et al., 2013).However, these changes are very difficult to predict: annual rainfall is predicted to increase in the Canadian Arctic (Ford et al., 2018), but this may be offset by increases in evapotranspiration and decreases in permafrost, leading to soil drying (Meredith et al., 2022).Another source of uncertainty is the extent to which the modelling assumptions of habitat suitability models reflect the reality.
For example, species-habitat relationships may change over time.Like all correlational habitat suitability models, our methods assume these relationships will be stable over time, but the relationship between species and their environment may change under scenarios of climate change.Our models also assume that these species are at equilibrium with their environment.Shorebird populations have declined and therefore some potentially suitable habitat may be unoccupied, which may result in habitat suitability models that underestimate the area of potential suitable habitat.The data for our analysis was collected over a 25 year time span during which shorebird distributions appear to be shifting in response to climate change (Anderson, Fahrig, Rausch, Martin, et al., 2023), increasing the uncertainty in the present day habitat suitability models upon which the future projections are built.
Despite the uncertainty, the consensus picture (Figure 5) suggests a northward and westward shift of breeding habitat for shorebirds Future projections of habitat suitability and species distributions are frequently proposed as something to consider in the design of climate resilient protected areas networks (Araújo, 2009;Ochoa-Ochoa et al., 2016).Countries around the world are pursuing protection of 30% of their land mass under the Global Biodiversity Framework (Xu et al., 2021).However, comprehensive information about the degree of uncertainty in these projections is needed in order to understand whether this is desirable or even possible.The probabilistic methods we used here are important for presenting a more honest representation of the uncertainty around future projections of habitat suitability and could inform the development of climate resilient protected areas networks.
Unexpectedly, we identified that the modelling decision generating the most uncertainty in our predictions was how candidate variables were selected to be included in the model, a source of uncertainty that has not been recognized in previous assessments of uncertainty in predictions of habitat suitability under climate change.
Our study also sheds new light on the potential impacts of climate change on Arctic-breeding shorebird species in Canada.While previous studies have described a significant risk of losing more than half of their suitable breeding habitat, our findings indicate that this prediction is accompanied by a high level of uncertainty.Explicitly accounting for uncertainty when projecting how habitat of vulnerable species like shorebirds will shift in response to climate change will be key for making conservation decisions that succeed across a range of potential futures.

ACK N OWLED G EM ENTS
We thank the field crews who conducted shorebird surveys for the
used SRES carbon emissions scenarios B2b and A2a; we used close modern equivalents, CMIP6 Shared Socio-economic Pathways SSP 245 ("middle of the road" emissions) and SSP 585 ("fossil-fueled development") emissions (O'Neill et al., 2017).The future projections for land cover included three possible scenarios for how the tree line could disperse northwards: a maximum shift of 5 km over 50 years, a maximum shift of 20 km over 50 years, or an equilibrium scenario where there is no maximum distance imposed on the movement of the tree line.
E 2 Mean probability distributions for proportional change in suitable habitat (future area/present area) from 216 configurations of 6 modelling decisions (pool of candidate variables, model algorithm, selection of variables, global circulation model, carbon emissions scenario, and tree line dispersal).Probability distributions were calculated based on 500 bootstrapped resamples of shorebird occupancy data from the Canadian Arctic, 1994-2018 (see Figure 1).The black line indicates no change in the area of suitable habitat between 2010 and 2075.Red dashed lines indicate half or twice as much suitable habitat in 2075 compared to 2010.Grey lines indicate the minimum and maximum estimates for the change in suitable habitat for each species.

F
I G U R E 3The proportion of the sum of squares attributed to different sources of uncertainty in ANOVA models of the proportional change in suitable habitat (a, b), the probability of predicting half the amount of suitable habitat (c), and the probability of predicting twice the amount of suitable habitat (d) in 2075 compared to 2010.In panel a, we compare the sum of squares attributed to model configuration, a factor with all 216 possible configurations of 6 modelling decisions, to the sum of squares attributed to species and the residual unexplained variation.In panels b-d, we break down the sum of squares attributed to model configuration to look at the uncertainty attributed to each of the 6 individual modelling decisions, with species as a random effect.
sources of uncertainty.Galbraith et al. (2014) used a qualitative risk assessment framework to evaluate the vulnerability of North American shorebirds to climate change.They predicted that 9 of 10 shorebird species from our study will lose at least 50% of their breeding habitat, and claimed high subjective confidence (>70% probability) in that prediction.Galbraith et al. (2014) may have higher confidence in predicting major losses of breeding habitat because their methods were based on expert opinion instead of modelling.

F
I G U R E 4 Across 10 Arctic-breeding shorebird species, the effects of six modelling decisions on the predicted probability that a species' area of suitable habitat in 2075 will be half or twice as much as in 2010.The box and whiskers show the median, 25th and 75th percentiles, and 1.5 times the interquartile range, with individual dots for values outside this range.The box plot for each modelling decision includes all other possible combinations of the other 5 modelling decisions.variable selection has not been previously identified as an important source of uncertainty in projections of habitat change.Dormann et al. (2008) concluded that the choice of variable selection was less important than the choice of modelling algorithm as a source of uncertainty in climate change projections, but their study only in the Canadian Arctic.Several areas that are currently considered good shorebird breeding habitat, including Victoria Island (Kitlineq), Southampton Island (Shugliaq), Prince Charles Island, and the Dewey Soper (Isulijarnik) Migratory Bird Sanctuary(Latour, 2008), were predicted to decline in the number of shorebird species they support.The areas where suitability for shorebirds will increase include Prince of Wales Island, and the Queen Elizabeth Islands, including Qausuittuq National Park and Quttinirpaaq National Park.Despite the risks of losing shorebird habitat, the northward shift we predict in suitable shorebird habitat suggests that Canada will be an important refuge for Arctic breeding shorebirds, because in many other regions of the Arctic there is no land farther north for these species to shift northwards into.
PRISM program, and the Polar Continental Shelf program for providing logistical support.We also thank the anonymous reviewers for their comments on the manuscript.CMA is supported through scholarships from the Natural Sciences and Engineering Research Council, Ontario Graduate Scholarships, and The W. Garfield Weston Foundation.PAS is supported through grants from the Natural Sciences and Engineering Research Council and ArcticNet NCE.