Rapid shifts in Arctic tundra species' distributions and inter‐specific range overlap under future climate change

The Arctic is one of the planet's most rapidly warming regions, with trends expected to intensify in the future. Projections of shifts in species distributional ranges under future climate change are thus far lacking for most vertebrate species using the Arctic tundra. Our aim was to assess possible climate‐induced changes in distributional ranges and inter‐specific overlap of an Arctic species assemblage within the world's largest land‐based protected area.


| INTRODUC TI ON
Climate is a major determinant of species distribution ranges (Louthan et al., 2015). Future changes in climate are expected to shift distribution ranges poleward (higher latitude) and/or upward (higher elevation) (Lenoir et al., 2008;Parmesan & Yohe, 2003).
Yet, due to regional variation in the magnitude and direction of climatic changes, impacts will vary between ecosystems (IPCC, 2013). Moreover, differences in physiological traits (e.g. thermal tolerance) between species can lead to divergent responses to global warming (Bennett et al., 2019). In agreement, species distributions have been shown to shift, contract or expand (Chen et al., 2011;Scheffers et al., 2016), which can lead to a restructuring of biotic communities over large areas (Lavergne et al., 2010).
The impacts of climate change are believed to be particularly pronounced for cold-adapted species residing in the Arctic, an area disproportionally affected by climate change relative to most other regions in the world Post et al., 2009). While the ecology and evolution of some resident Arctic species are already being affected by a warming climate Taylor et al., 2020), most of the evidence is based on marine (Clairbaux et al., 2019;Kovacs et al., 2011;Neukermans et al., 2018) and invertebrate Nielsen & Wall, 2013) species. Projections of distributional range changes by an assemblage of vertebrate species found on the Arctic tundra and possible implications for interspecific overlap are currently lacking. Such knowledge is crucial given that the Arctic tundra has low species richness with relatively few species interactions, where even small changes in inter-specific range overlap can have important implications for terrestrial foodweb dynamics (Schmidt et al., 2017). Estimating climate-induced shifts in distributions of Arctic species is therefore urgently needed to facilitate assessments of potential implications for ecosystem functioning (Pecl et al., 2017) and to inform possible conservation strategies (Elsen et al., 2020).
Following the general notion that climate change impacts are more pronounced in areas with higher levels of warming, such as the Arctic (Chen et al., 2011;IPCC, 2013), we expected species distribution ranges to expand northwards in latitude and up the elevation gradient over time (P 1 ). We also expected the speed of northward and upward shift to depend on the severity of climate change (i.e. RCP scenarios) (P 2 ). Moreover, we expected species-specific responses to climate change to alter inter-specific range overlaps in space and time (P 3 ).

| Study area and location data
Data for this study were collected in the Northeast Greenland National Park (Figure 1). The size of the park is 972,000 km 2 including the interior ice sheet (ca 70% of the area). Apart from a few permanent, small meteorological and military outposts, the national park has no permanent human residents. Only researchers and few tourists visit the park during summer months, and anthropogenic disturbance is thus considered low. Large areas of the ice-free part of the national park are patrolled by the Sirius Dog Sled Patrol, a special unit within the Danish defence. The patrol consists of 6 dog sled teams, each conducting one reconnaissance trip in spring and one in autumn of each year. The exact routes are classified, but trips cover the same geographical areas, each covering over a thousand kilometres across several months. During reconnaissance trips, vertebrate wildlife encounters are registered, specifying date, species, location and group size.
Here, we used species observation data collected during 1979-2013. Observations collected during 2001-2013 were available in a digital format while observations from 1979 to 2000 were available on paper only. We first digitized these observations, considering only records with known species identity, date and location. The locations of most observations were provided in lat/long, but in some cases, only local place names were given, which were subsequently converted into lat/long using GIS (https://data.geus.dk/geusm apmor e/stedn avned b/lista ll.jsp). We could not quantitatively assess the accuracy of the location data, but we assume location error of most observations to be <20 km, which is smaller than the resolution of the bioclimatic variables used in the analyses (see next section). We disregarded information of group size, as this was inconsistently registered for most species.
For our analyses, we focussed on location data of eight vertebrate species that were easily identified by Sirius Dog Sled Patrol personnel due to their size, conservation status or charismatic na- as lemming, were not considered for the analyses as they were less likely to be observed and recorded, even though they may be important to tundra ecosystem functioning . The location data reflect species presences only, as there was no data recorded for true absences of species over time and space. We therefore also considered observations of animal tracks (<1% of location data) in the analyses. Part of the muskox location data collected by the Sirius Dog Sled Patrol has been used previously to assess biotic and abiotic effects on temporal population trends (Forchhammer & Boertmann, 1993).

| Species distribution modelling
Species distribution ranges were estimated using the maximum entropy machine learning algorithm (MaxEnt; Phillips et al., 2006). MaxEnt relates location data with environmental or climatic background data to produce spatially explicit predictions of species' occurrence probability . Although a wide range F I G U R E 1 Map of Greenland with the Northeast Greenland National Park outlined in green and the Greenland ice sheet in blue. Also shown are the study species: Arctic fox (Vulpes lagopus), Arctic hare (Lepus arcticus), Arctic wolf (Canis lupus arctos), muskox (Ovibos moschatus), polar bear (Ursus maritimus), rock ptarmigan (Lagopus muta), snow bunting (Plectrophenax nivalis) and snowy owl (Bubo scandiacus) and their respective sample sizes (N is the number of observations over the period 1979-2013) as used in the species distribution models. Arrows indicate the most common, but not exclusive, interactions between carnivorous and herbivorous species. Note that arrows and artwork are for descriptive purposes only and are thus not to scale of species distribution modelling (SDM) techniques exists, which can be combined into a single ensemble modelling framework (Forester et al., 2013;Naimi & Araújo, 2016), we opted for MaxEnt as the preferred modelling approach as it is particularly suited for presenceonly data with relatively small sample sizes (Elith et al., 2006;Phillips et al., 2017;Wisz et al., 2008), which was the case in our study ( Figure 1, Figure S1.1 in Appendix S1). We fitted separate models for each species using the procedure outlined below. For more details, we refer to the Overview, Data, Model, Assessment and Prediction (ODMAP) protocol (sensu Zurell et al., 2020) on model development, testing and evaluation in Appendix S1.
Presence data in the MaxEnt models were species' locations collected in the study area during 1979-2013 ( Figure S1.2 in Appendix S1).
Background points (10,000) were randomly sampled for each species from within the vegetated areas in the Northeast Greenland National Park ( Figure S1.3 in Appendix S1), effectively removing all areas with land ice. Moreover, we constructed spatial sampling bias files for each species separately by computing Gaussian kernel density rasters of all sampling locations (Brown et al., 2017). Sampling bias files ( Figure   S1.4 in Appendix S1) were subsequently included in the speciesspecific MaxEnt models to up-weight presence-only data points with fewer neighbours in the geographic landscape and to restrict background points to geographic areas where species occurrences were found (Merow et al., 2013;Phillips et al., 2009).
Environmental data considered were bioclimatic variables (n = 19) freely available through the CHELSA database (https://chels a-clima te.org/) at a 1 km resolution (Karger et al., 2017). The bioclimatic raster data were downloaded for the periods 1979-2013 (current period) and 2061-2080 (future period). For the latter, we considered data produced by three global circulation models (GCMs) using three scenarios of representative concentration pathways (RCPs). Specifically, the three selected GCMs are part of the CMIP5 collection of model runs used in IPCC's 5th Assessment Report (IPCC, 2013) and were CM5A-LR (Persechino et al., 2013), ESM-MR (Giorgetta et al., 2013) and MIROC5 (Watanabe et al., 2010). Each GCM provides data for scenarios of RCP 2.6, 4.5 and 8.5 (see Appendix S1 for a description of each RCP scenario). All bioclimatic raster data were upscaled to 20 km resolution to minimize possible effects of location inaccuracy and clipped to the vegetated areas in Northeast Greenland National Park ( Figure S1.3 in Appendix S1).
Because multicollinearity among bioclimatic predictors was high (variance inflation factor>100 for some variables), and we had no a priori knowledge of which predictor variables were most influential in explaining species-specific probability of occurrence, we employed a completely data-driven variable selection approach through the "SDMtune" package in R (Vignali et al., 2020). This approach selects the variable that best fits the data (using Akaike information criterion value corrected for small sample sizes: AICc) among those that are highly correlated. Starting from a trained model with all variables included, the function checks if the variable ranked as the most important (i.e. highest per cent contribution) is correlated with any of the other variables, here using a threshold of Spearman's rho ≥0.7.
If this was the case, a leave-one-out jackknife test was performed, starting with the full model, and among all correlated variables, the one variable that decreased model fit on the training data set the most was discarded. A new model was then trained without this variable, which was again checked for highly correlated variables.
To protect against overfitting and to reduce model complexity, MaxEnt uses regularization multipliers (RM) (Phillips et al., 2006).

RMs give a penalty for each term included in the model and for
higher weights given to a term. Here, we tested different settings of RM using the range of 0.5-5.0 in increments of 0.5 for each feature class through the "ENMeval" package in R (Muscarella et al., 2014).
Moreover, we restricted all features to "linear" and "quadratic" functions to avoid overly complex response curves that would be hard to explain ecologically. For each species, the optimal model was selected based on goodness of fit by identifying the candidate model with the lowest AICc. Predictive performance of each speciesspecific optimal model was evaluated by calculating the area under the curve (AUC) of the receiver operating characteristic (ROC) value.
AUC values typically range from 0.5 for a model that shows no better than random discriminatory ability, to a theoretical value of 1 for perfect discriminatory ability (Elith et al., 2006). Autocorrelation of location data was considered low (Figure S1.6 in Appendix S1).
Species-specific distribution maps were created by stacking the bioclimatic raster layers into a multi-layered raster and predicting, from the optimal MaxEnt models, the probability of occurrence in

| Shifts in distribution ranges and interspecific overlap
In order to investigate how changes in bioclimatic conditions might shift distributional ranges, we contrasted the predicted probability of occurrence, as derived from the complimentary log-log (cloglog) output produced by the optimal MaxEnt model, between the current and future periods. To do so, we considered three complementary SDM thresholds (Liu et al., 2013) including Kappa (the value of the probability of occurrence at which Kappa is highest), MSSS (the value of the probability of occurrence at which the sum of the sensitivity (true positive rate) and specificity (true negative rate) is maximized) and P10 (the value of the probability of occurrence for the lowest 10% of occurrence records). In general, the Kappa threshold was most restrictive as it identified areas with relatively high probability of occurrence (i.e. core distribution area). The MSSS threshold identified areas above a moderate probability of occurrence, while the P10 threshold included most areas above a relatively low probability of occurrence across species distributions (Table S1.1 in Appendix S1).
Directional shifts under recent and future climatic conditions (considering the three RCP scenarios) were then assessed quantitatively by calculating the difference in latitude (using the geometric mean, with 1° corresponding to 111.6 km) and elevation (mean m.a.s.l.) between periods (a difference of five decades) and for each SDM threshold. Elevation was extracted from a digital elevation model (20 km resolution) developed by the Shuttle Radar Topography Mission, which is freely available through the WorldClim2 database (http://www.world clim.com/). Moreover, we computed the absolute change in total area size (km 2 ) of species ranges and the level of clustering (unitless) (Clark & Evans, 1954). For the latter, we followed For each range shift metric, we used a paired sample t test to detect significant (p < .05) differences in mean species distributions between the current and future period (normality of data confirmed with Shapiro-Wilk test: p > .05) with a separate test for each RCP and SDM threshold considered. Finally, we calculated the proportion of inter-specific overlap, and the changes therein between recent and future species distributions, based on Schoener's D niche overlap index (Schoener, 1968), again separating each RCP and SDM threshold considered.

| Predictive performance and variable importance
Predictive performance of the species-specific MaxEnt models was considered high with a mean AUC >0.75 for all optimal models (Table   S2.1 in Appendix S2). The importance of bioclimatic variables in the models and their influence on the probability of occurrence differed between species (see Figure S2.1 for variable importance and Figure   S2.2 for response curves). The most consistent and important predictor influencing the probability of occurrence across all species was mean temperature (°C) of the warmest quarter (summer) with the probability of occurrence typically increasing with increasing temperature. These forecast maps suggest that the probability of species occurrences will increase and expand geographically and particularly with severity of RCP scenario, a pattern that applied to all species. The maps also suggest that some GCMs produced more spatially restrictive distribution ranges under the same RCP scenario than did others. To accommodate this variance in further analyses, we calculated the mean and standard deviation of the probability of species occurrence for each RCP by aggregating data from the GCM-and speciesspecific MaxEnt models ( Figures S2.11-S2.12 in Appendix S2).    (Table S2.2b in Appendix S2). The total area size of projected core ranges more than doubled over time for nearly all species and RCP scenarios ( Figure 4c, Table S2.2c in Appendix S2), with an average increase of 18,700 km 2 /decade across species under RCP 8.5 (t = 7.47, df = 7, p < .001). While we found a general tendency for predicted distribution ranges to become less clustered (i.e. more dispersed) compared with the current situation (Figure 4d), a statistical difference was only found for the MSSS threshold under the scenario RCP 8.5 and for the P10 threshold across RCP scenarios (Table S2.2d in Appendix S2). As such, predicted core ranges were found to largely maintain their current degree of clustering (Figure 4d, Table S2.2d in Appendix S2).

| Changes in inter-specific range overlap
The projected shifts in species-specific distributions over time led to distinct changes in inter-specific range overlap, though patterns varied between the RCP scenarios and thresholds considered ( Figure 5).

F I G U R E 2
Maps of species-specific probability of occurrence for the period 1979-2013 based on the optimal MaxEnt models using location data collected in the Northeast Greenland National Park during the same period. Predicted values are the cloglog output of the species-specific MaxEnt model with values ranging from 0 to 1 depicted by a blue to green scale. Note that values are only provided for vegetated areas in the park, excluding land ice and areas with novel bioclimatic conditions as identified through species-specific multivariate environmental similarity surfaces analyses using data from 1979-2013 and 2061-2080 In general, proportional range overlap between most species' combinations increased with severity of the climate change scenarios.
However, the increase was most pronounced for the threshold P10.
Proportional range overlap for the thresholds Kappa and MSSS was projected to decrease or remain stable between most species' combinations under scenario RCP 2.6 and 4.5 and increase under RCP 8.5 only ( Figure 5).

| D ISCUSS I ON
This study provides a comprehensive overview of how projected climate change scenarios for north-east Greenland may shift species' distributions and subsequently inter-specific overlap of eight Arctic vertebrate species within the world's largest land-based national park. Following our expectation (P 1 ), we found strong indications for northwards range shifts by all species when contrasting probability of occurrences between recent  and future (2061-2080) bioclimatic conditions. Similarly, projected distribution ranges were found to move up the elevational gradient (also supporting P 1 ).
Both of these directional shifts were apparent across SDM thresholds considered, and clearly accelerated with the severity of the RCP scenarios (supporting P 2 ). In fact, the range shifts across RCP scenarios found here were all greater than reported in a previous meta-analysis (16.9 km/decade over the latitudinal gradient) that was based on a wider array of taxa and ecosystems from across the globe (Chen et al., 2011). As such, our Arctic case study underlines the notion that climate change impacts are most pronounced in areas with higher levels of warming, leading to accelerated shifts in species' ranges towards the poles (Chen et al., 2011;IPCC, 2013). Importantly, we also found a clear expansion of species' distribution ranges under future conditions, with core ranges more than doubling in size for nearly all species and RCP scenarios. Together with a general tendency for less clustered distribution ranges under the RCP scenarios considered, we conclude that future climatic conditions in vast areas in north Greenland may not have such a destressing impact on the distribution of cold-adapted species than perhaps previously suspected, at least not within this century.
Species-specific responses to changing climatic conditions are often detected and likely driven by differences in physiological traits (e.g. thermal tolerance), which can culminate into a restructuring of biotic communities over large areas (Lavergne et al., 2010). In agreement, the importance of bioclimatic variables included in our MaxEnt models and their influence on the probability of occurrence differed between species, driving speciesspecific changes between current and projected distributions. The Results are provided for three SDM thresholds: Kappa, MSSS and P10. Significant differences between current and future values, as determined by paired sample t tests, are indicated with brackets and stars where * indicates p < .05, ** indicates p < .01 and *** indicates p < .001 (see Table S2.2 in Appendix S2 for full test outputs). Nearest neighbour index values closer to 1 indicate reduced clustering and more evenly dispersion of occurrence. Values were derived based on the species-specific optimal maximum entropy (the complimentary loglog output) models as identified by the Akaike information criterion value corrected for small sample sizes 2.6). In contrast, we found an increase in overlap between all species combinations under the most severe scenario that would limit global temperature rise to 4°C in year 2100 (RCP 8.5). Overall, our findings do not lend much support to the notion that that future climate change in this region will lead to a complete rearrangement of trophic interactions and food-web dynamics, as hypothesized elsewhere (Legagneux et al., 2014), though further studies are required to fully address this issue.
Efforts to predict shifts in species distribution ranges and overlap are critical to inform management and conservation initiatives (Elsen et al., 2020), but estimating climate change effects with SDMs alone remains challenging (White et al., 2018). Our correlative analyses focus purely on estimating possible shifts in distribution ranges based on gradual changes in climatic conditions and ignore possible biotic interactions. Indeed, our approach follows the Eltonian noise hypothesis (ENH), which proposes that although biotic interactions   (Stewart et al., 2018). Doing so at the landscape scale, as would be required in our forecasting analyses, remains a complicated task . This is due to the multitude of biotic and abiotic interactions that vary over space, as well as time lag effects involved between loss of snow/ice cover and the build-up of soil depth and quality that is needed for vegetation growth (Jobbágy & Jackson, 2000). Using bioclimatic variables as surrogates for actual limiting resources, such as forage availability and quality, is often the only reliable alternative to assess possible changes in species ranges under future global warming. We thus interpret our results under the assumption that biotic factors can correlate closely with abiotic variables and that our SDMs capture an important part of the biotic processes involved in species distributions.
Systematically collecting long-term location data of wildlife across large areas in the high-Arctic is extremely challenging and expensive, and thus rare. The here analysed data set is likely the most extensive that is currently available for within Greenland, possibly even within the Arctic region. Nonetheless, some challenges in the data set required consideration so as to reduce prediction uncertainty, which is often neglected in large-scale SDM studies investigating climate change impacts (Beale & Lennon, 2012). First, location accuracy of the presence data was unknown, yet likely poor in comparison with, for example GPS tracking or camera trap data. We attempted to minimize this uncertainty related to location error by using coarse grid cells of bioclimatic raster data (20 km 2 ).
Next, an important assumption of SDM studies is that sampling of location data is adequate and representative, which we could not quantify here given the rather opportunistic nature of how species observations were collected using unknown routes of the Sirius Dog Sled Patrol. To account for this issue as much as possible, we incorporated spatial sampling bias files in the species-specific MaxEnt models, which is an established method to restrict background points to geographic areas where species occurrences were found, leading to more realistic predictions (Phillips et al., 2009). We also tailored the entire analytical procedure to increase reliability of model output as much as possible by, for example incorporating a wide range of emission scenarios, excluding areas with novel bioclimatic conditions in the future and limiting overparameterization through extensive MaxEnt model pruning (Santini et al., 2021).
We also explicitly quantified a known source of uncertainty in climate change studies, namely the GCMs considered in forecasting (Thuiller et al., 2019). Finally, it is important to reiterate that our study system is a remote and protected area, currently unaffected by human disturbance. Anthropogenic influences on contemporary species distributions can therefore be considered low, which eliminates an otherwise known bias in forecasting climate change impacts on species ranges (Faurby & Araújo, 2018). However, as the Arctic climate warms and the ice continues to melt, human activity in the region (e.g. oil and mineral extraction) is steadily growing and will likely play an increasingly important role in shaping future distribution ranges of wildlife species resident to the Arctic (Johnson et al., 2005). Despite these caveats, and based on a purely bioclimatic assessment of species-specific shifts in distribution and overlap under future climate change scenarios, our results clearly show that ongoing climate warming is likely to have a strong impact in the Arctic, with rapid and directional shifts in species' distribution ranges towards higher latitudes and elevations. To what extent associated changes in inter-specific overlap of distribution ranges under future conditions will impact local food-web dynamics and ultimately ecosystem functioning remains to be studied in more detail. Our results therefore emphasize that understanding climate change effects on terrestrial food-web dynamics in the Arctic should remain a research priority.

ACK N OWLED G EM ENTS
We thank all sled patrol teams whose dedicated efforts over the years have made this study possible. We also thank the Joint Arctic Command for making the data available to us. We also acknowledge the Editor and an anonymous reviewer for their constructive comments on a previous draft, which greatly improved the paper.